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A Survey of the Core-Congruential Formulation for 
Geometrically Nonlinear TL Finite Elements 

C. A. FELIPPA*, L. A. CRIVELLI ** and B. HAUGEN *** 


SUMMARY 

This article presents a survey of the Core-Congruential Formulation (CCF) for geometrically 
nonlinear mechanical finite elements based on the Total Lagrangian (TL) kinematic description. 
Although the key ideas behind the CCF can be traced back to Rajasekaran and Murray in 1973, it 
has not subsequently received serious attention. The CCF is distinguished by a two-phase devel- 
opment of the finite element stiffness equations. The initial phase develop equations for individual 
particles. These equations are expressed in terms of displacement gradients as degrees of freedom. 
The second phase involves congruential-type transformations that eventually binds the element 
particles of an individual element in terms of its node-displacement degrees of freedom. Two ver- 
sions of the CCF, labeled Direct and Generalized, are distinguished. The Direct CCF (DCCF) is 
first described in general form and then applied to the derivation of geometrically nonlinear bar, 
and plane stress elements using the Green- Lagrange strain measure. The more complex Gener- 
alized CCF (GCCF) is described and applied to the derivation of 2D and 3D Timoshenko beam 
elements. Several advantages of the CCF, notably the physically clean separation of material and 
geometric stiffnesses, and its independence with respect to the ultimate choice of shape functions 
and element degrees of freedom, are noted. Application examples involving very large motions 
solved with the 3D beam element display the range of applicability of this formulation, which 
transcends the kinematic limitations commonly attributed to the TL description. 


1. INTRODUCTION 

There is an elegant Total Lagrangian (TL) formulation of geometrically nonlinear mechani- 
cal finite elements that has received little attention in the literature. This will be referred to 
as the Core-Congruential Formulation , or CCF, in the sequel. The key concepts, presented 
by Rajasekaran and Murray 1 in 1973, evolved from the analysis and reinterpretation of the 
pioneer work of Mallet and Marcal, 2 as well as Murray’s previous work in geometrically 
nonlinear finite element analysis. 3 The discussion of Reference 1 by Felippa 4 provided para- 
metric expressions for the stiffness matrices that appear at various levels of the discrete 
governing equations. This work originated what is called here the Direct Core Congruential 

Formulation, or DCCF. 

In 1987 a course in nonlinear finite elements offered by the first author presented the 
derivation of several elements using the DCCF. Preparation of homework assignments and 

* Department of Aerospace Engineering Sciences and Center for Aerospace Structures, Uni- 
versity of Colorado, Boulder, CO 80309-0429, USA 

** Hibbit, Karlsson k Sorensen Associates, 1080 Main St., Pawtucket, RI 02860, USA 

*•* Division of Structural Engineering, The Norwegian Institute of Technology, N-7034 Trond- 
heim, Norway 


1 


feedback from students in this and follow-up offerings helped to streamline the material. 
Subsequently Crivelli s doctoral thesis 2 * * * 6 used the CCF in the systematic development of a 
three-dimensional nonlinear Timoshenko beam element capable of undergoing arbitrarily 
large rotations. Challenges posed by this application pushed this formulation beyond 
frontiers hitherto deemed impassable by a TL element with rotational degrees of freedom. 
This development was summarily reported in a survey article by Felippa and Crivelli 6 and 
explained in more detail in a subsequent paper by Crivelli and Felippa. 7 * 

A lesson gained from this research is that, when dealing with 3D finite rotations, the CCF 
should be applied in a staged fashion that allows the systematic examination of additional 
terms arising in the transformations to physical degrees of freedom. That transformation 
methodology gave rise to what is here called the Generalized CCF, or GCCF. 

Both DCCF and GCCF share the same u divide and conquer” philosophy. However, the 
core equations as well as subsequent steps that transform those equations to physical 
freedoms vary in complexity. To simplify the exposition while focusing on the essential 
aspects, Sections 3 through 7 focus on the DCCF. Examples of application to elements 
amenable to the direct treatment are presented. The GCCF is discussed in Sections 8 
through 10, and illustrated with applications to 2D and 3D beam elements. 

REMARK 1.1. Several authors have expressed the belief that the approximation performance of 
TL- based elements degrades beyond moderate rotations, and an updated Lagrangian or corota- 
tional description is necessary for handling truly large motions. For example, in 1986 Mathiasson, 
Bengtsson and Samuelsson concluded that “The TL formulation can only be used in problems 
with small or moderate displacements.” More recently Bergan and Mathisen 9 voice a similar 
opinion: “it is commonly known that in a step by step TL formulation artificial strains easily 
arise in beam elements due to non homogeneities in the displacement expansions in transverse 
and longitudinal directions.” Our experience shows that such limitations are not inherent in the 
TL description but instead emerge when a priori kinematic approximations are made to simplify 
element derivations. The 3D beam element just cited exhibits computational and approximation 
performance for very large rotations comparable to those based on the co-rotational and Updated 
Lagrangian descriptions while retaining certain advantages listed in the Conclusions. 

2. OVERVIEW 


2.1 Basic Concepts 

The original development of the CCF was concerned with the construction of TL stiffness 

matrices for geometrically nonlinear analysis through the congruential-transformation pat- 

tern 

Rlevel = f G T S level G ^ ( 2 . 1 ) 

J Vo 

where S is the core stiffness matrix, K the physical stiffness in terms of the nodal degrees of 
freedom v, G a core-to-physical-freedom transformation matrix assumed to be independent 
of v, V 0 the appropriate reference integration volume, and in which “level” identifies the 
governing equation level at which the stiffness matrix is used. 
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The three variational levels of interest in practice are: energy (level 0), force equilibrium 
(level 1), and first-order incremental equilibrium (level 2). Qualifiers “residual-force” and 
“secant-stiffness” are also used for level 1, and “tangent-stiffness” used for level 2. 

The core stiffness matrix is expressed in terms of the displacement gradients at each ma- 
terial point. Displacement gradients g make a better choice of core variables than finite 
strains because for elements with translational degrees of freedom (DOFs) they can be 
expressed linearly in terms of node displacements v as g = Gv, a property that validates 
(2.1) for all levels. As discussed below, such elements fall under the purview of the Direct 
CCF. 

The qualifier “core” emphasizes the goal of independence of S level with respect to dis- 
cretization decisions such as element geometry, shape functions, and choice of nodal de- 
grees of freedom. Such a dependence is introduced by the congruential transformation 
indicated in (2.1) and the integration over the element volume. 

2.2 Direct and Generalized CCF 

The basic schematics of the CCF, mathematically expressed through (2.1), may be dia- 
grammed as 


Core 


Congruential 


Physical-DOF 

Stiffness 


Transformation 


Stiffness 

Equations 


Equations 


Equations 


But this panoramic view needs to be rendered more precise. If the relation between core 
DOFs (the displacement gradients g) and the physical DOFs (the node displacements v 
of a finite element model) is linear , these transformations do not depend on level: 


(0) Core Energy Stiffness 

(1) Core Secant Stiffness 
(1 ) Core Internal Force 

(2) Core Tangent Stiffness 


Congruential 

Transformation 

Equations 


Physical Energy Stiffness 
Physical Secant Stiffness 
Physical Internal Force 
Physical Tangent Stiffness 


In this diagram, numbers annotated within the “core box” denote the variational level 
of the governing equation in use. Internal force and secant stiffness are two alternative 
governing-equation expressions at level 1. The energy level (level 0) may also be expressed 
in several ways, but this is not shown in the diagrams to reduce clutter. Under the afore- 
mentioned assumption we obtain the Direct Core Congruential Formulation, or DCCF. 

If the relation between displacement gradients g and node displacements v is nonlinear, 
the transformations sketched above are not only more complex but depend on variational 
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level and possibly the expression form used within a level. This complication arises when 
elements with rotational degrees of freedom such as beams, plates and shells are considered. 
It gives rise to the Generalized Core Congruential Formulation, or GCCF. 

Two variants of the GCCF may be distinguished. If the relation between g and v is 
nonlinear but algebraic , the transformation equations do vary with level but in principle 
are still possible as illustrated in the following diagram. 


( 0 ) Core Energy Stiffness 

( 1 ) Core Secant Stiffness 


S- Congruential 
Transformation 
Equations 

=> 

Physical Energy Stiffness 
Physical Secant Stiffness 

( 1 ) Core Internal Force 
( 2 ) Core Tangent Stiffness 

=> 

T-Congruential 

Transformation 

Equations 

- 

Physical Internal Force 
Physical Tangent Stiffness 


Here “T- Congruential” and “S- Congruential” are abbreviations for “Tangent Congruen- 
tial” and “Secant-Congruential,” respectively. Such a distinction is elaborated upon in 
Section 8. 

If the relation between g and v is nonlinear and can be expressed only in non-integrable 
differential form, the “Secant Transformation Equations” of the preceding diagram do not 
generally exist, and the diagram must be truncated: 


( 0 ) Core Energy Stiffness 

( 1 ) Core Secant Stiffness 


( 1 ) Core Internal Force 

( 2 ) Core Tangent Stiffness 


T-Congruential 

Transformation 

Equations 


Physical Internal Force 
Physical Tangent Stiffness 


These two variants of the GCCF are called Algebraic GCCF and Differential GCCF and 
denoted by acronyms AGCCF and DGCCF, respectively, in the sequel. The main dis- 
tinction between AGCCF and DGCCF is that it makes no sense to talk about missing 
quantities, such as the physical secant stiffness, with the latter. 

The original development of the CCF outlined in the Introduction focused on elements 
with translational-degree-of-freedom configurations. For such elements the Direct form 
of the CCF, or DCCF, is sufficient. Sections 3 through focus on that form, leaving the 
development and application of the GCCF to Section 8 and following ones. 
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2.3 The CCF Philosophy: Divide and Conquer 


The CCF derivation of the finite element equations naturally reflects the outlined frame- 
work. It proceeds through two phases: a core phase followed by a transformation phase. 
In the initial phase core energy, secant and tangent stiffness matrices as well as internal 
force vectors are obtained. These matrices and vectors pertain to individual particles. For 
the stiffness matrices they are collectively represented by the term S levet in (2.1). 

The key goal is to try to make such core equations as independent as possible with re- 
spect to finite-element discretization decisions such as element geometry, shape functions, 
selection of nodal degrees of freedom and (in the case of rotational DOFs) rotational 
parametrizations. To emphasize this independence, the term core was coined. Complete 
independence is in fact achievable if the relation between displacement gradients g and v 
is linear, which characterizes the DCCF. The goal has to be tempered if the relation is 
nonlinear because dependencies may arise at the tangent stiffness level. Such dependencies 
create the so-called complementary geometric stiffness terms, which are characteristic of 
elements that fall under purview of the GCCF. 

In the transformation phase, these core forms are transformed to physical DOFs, i.e. 
element node displacements. The transformation may be done directly for simple elements 
and in multistage fashion for complex ones. In particular, multistage transformations 
are recommended for elements that require the Differential GCCF such as 3D beam and 
shell elements. In this case the transformation phase is decomposed into transformation 
stages that progressively “bind” particles into lines, areas or volumes through kinematic 
constraints, and eventually link the element domain to the nodal degrees of freedom. 
Decisions such as the choice of specific parametrizations for finite rotations may be deferred 
to final stages. 

What are the differences between the CCF and the more conventional Total Lagrangian 
formulation of nonlinear finite elements? If kinematic exactness is maintained throughout, 
the final discrete equations are identical. This is shown in Appendix 1 for the DCCF applied 
to continuum elements. But in geometrically nonlinear analysis approximations of various 
kinds are common, especially in structural elements with rotational degrees of freedom 
such as beams, plates and shells. In the conventional formulation it is quite difficult to 
assess a priori the effect of seemingly innocuous approximations “thrown into the pot,” 
and a posteriori exhaustive testing of complex situations becomes virtually impossible. 
Sample: how does the neglect of higher order terms in the axial deformation of a spinning 
3D beam affects torsional buckling? 

The staged approach recommended for the GCCF permits a better control over such 
assumptions. The core equations are physically transparent, clearly displaying the effect of 
material behavior, displacement gradients and prestresses. In the ensuing transformation 
sequence the origin of each term can be accurately traced, and on that basis informed 
decisions on retention or dropping made. This process can be aided by computer by 
testing subproblems that isolate the physics modeled by specific terms. 
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From this discussion it follows that, from the standpoint of element development, evalu- 
ation and testing, the most significant advantage that can be claimed for the CCF is the 
clean separation of physical effects. The importance of this factor should not be underes- 
timated, because physical transparency is the key to success in nonlinear analysis. 

3. HISTORICAL BACKGROUND 

In 1968 Mallet and Marcal 2 attempted to establish a standard nomenclature for geomet- 
rically nonlinear finite element structural analysis based on the Total Lagrangian (TL) 
kinematic description. Consider a discrete, finite element model of a static structural sys- 
tem under dead loading with nodal displacement degrees of freedom collected in array 
v. Displacements are measured from a fixed reference configuration Co to a current con- 
figuration C. The virtual- work conjugate forces, independent of v, are collected in array 
p. The system has a total potential energy function J = U - W that is the difference 
between the strain energy U and the loads potential W = p T v. The residual node forces 
are r = dJ/dv, and the symbol A denotes increment associated with the variation of the 
current configuration. (In keeping up with the spirit of Reference 2 actual variations are 
used in this Section rather than virtual ones; the latter are identified by the usual 6 prefix.) 

Mallet and Marcal expressed the total potential energy, the residual (force-balance) equi- 
librium equations, and the incremental equilibrium equations as follows: 


J = u- W = fv r [Ko + JN, + |N 2 ]v-p T v, 

(3.1) 

dJ 

r = = [K 0 + + |N 2 ] v - p = 0, 

(3.2) 

Ar = [Ko + Ni + N 2 ] Av — Ap = 0. 

(3.3) 


Here K 0 is the linear stiffness matrix evaluated at the reference configuration, whereas 
Nj and N 2 are nonlinear stiffness matrices , also evaluated at the reference configuration, 
that depend linearly and quadratically, respectively, on the node displacements v. The N 
matrices were said “to repeat” in the foregoing expressions. (This old notation has not 
survived; presently symbol N is most commonly used to identify matrices of element shape 
functions.) 

Five years later Rajasekaran and Murray 1 examined more critically the structure of the 
matrices that appear in the above equations. In that investigation they chose to start 
from the “core” stiffness matrices corresponding to K, Ni and N 2 expressed in terms 
of displacement gradients, and in doing so laid down the main idea of the CCF. Working 
with specific elements they showed that the nonlinear stiffness matrices Nj and N 2 are not 
uniquely determined. Indeed (3.1)-(3.3) as written are unique only for a single degree of 
freedom. They did not present, however, a general expression valid for arbitrary elements. 
This was partly done by Felippa, 4 who in the discussion of Reference 1 considered again 
those equations, rewritten here in a more general and compact form: 

J = Jv^v + (p° - p) T v, (3.4) 
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r = K r v + p° — p = f — p = 0, 
Ar = K Av — Ap = 0, 


(3.5) 

(3.6) 


in which the notation of this paper — rather than that of Reference 4 — is used. Here K 17 , 
K r and K denote the energy , secant and tangent stiffness matrices, respectively. (Energy 
and secant stiffnesses are not denoted by K e and K J because such symbols are used for 
other purposes in the finite element course noted in the Introduction.) In addition, p° is 
the prestress force vector , which vanishes if the reference configuration is stress free and 
weis omitted in that discussion, 4 and f = K r v-f p° is the internEil force vector. The tangent 
stiffness is of course fundamental in incremental-iterative solution methods and stability 
analysis, while the secant stiffness (by itself or in the internal- force form K r v + p° ) is 
important in pseudo-force methods. The energy stiffness enjoys limited application per se 
but has theoretical importance as source for the other two. 

In linear problems = K r = K = Ko and the three stiffness matrices coalesce. But 
in nonlinear problems not only do the matrices differ but, as shown in the next section, 
K 17 and K r may involve arbitrary scalar coefficients. Such parametrized expressions were 
given by Felippa 4 under the following restrictions: 

(Rl) K r is symmetric. 

(R2) The reference configuration is stress free. 

(R3) The finite strain measure is quadratic in the displacement gradients. 

(R4) The transformation between core and physical freedoms is linear. 

The following treatment eliminates restrictions (Rl) and (R2) altogether, and the other 
two selectively. It should be noted that restriction (R4) is the condition that, with present 
terminology, characterizes the DCCF. 

4. CORE STIFFNESS EQUATIONS 
4.1 TL Description of Particle Motion 

A conservative, geometrically nonlinear structure under dead loading is viewed as a con- 
tinuum undergoing finite displacements u. These displacements are measured from a fixed 
reference configuration C 0 to a variable current configuration C. No discretization into 
finite elements is implied at this stage. We confine our attention to the case in which the 
material behavior stays within the linear elastic range, thus implying small deformational 
strains but arbitrarily large rotations. Corresponding points or particles in the reference 
and current configuration are referred to a fixed Cartesian coordinate system and have 
the coordinates Xi and X{ ( i = 1, . . . rid), respectively, where nj is the number of space 
dimensions. The displacement field components are U{ = — Xi. 

Let the state of strain at a particle in the current configuration be characterized by n s 
strains e, (i = 1,2, ... n,) collected in an array e, and let the corresponding conjugate 
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stresses be .s, (i — 1,2,... n s ), collected in an array s. Using the summation convention 
the elastic stress-strain relations are written 

Si = s°i + Eijej, with E ig = Eji, or s = s° + Ee, (4.1) 

where s ® are stresses in the reference configuration (stresses that remain if e, = 0, also 
called piestresses) and Eij are elastic moduli arranged as a n 3 x n s square array in the 
usual manner. 

Let J, U, W, $ and T denote the analogues of J, U , W, p, f and r, respectively, at 
the particle level. (The first three acquire the meaning of energy densities, whereas is 
a dead-loading body force density independent of u.) The strain energy density can be 
expressed as 

U = e,s° + \e{Eijej = e T s° + |e T Ee. (4.2) 

The total strain energy U is obtained by integrating (4.2) over the structure volume: 
U fy 0 M dV~ the integration taking place — as can be expected in a TL description — 
over the reference configuration geometry. 

Next, introduce the n g displacement gradients g m n = du m /dX n . These are subsequently 
identified as g t (i = 1,2,... n g ) so they can be conveniently arranged in a one- dimensional 
array g. Following Rajasekaran and Murray 1 and Felippa 4 assume that the strains e; are 
linked to the displacement gradients through matrix relations of the form 

e, = hf g + |g T Hi g , * = 1, 2, . . . n a (4.3) 

where h, and H, are arrays of dimension n g x 1 and n g x n g , respectively, with H, sym- 
metric. In the original References 1 ’ 4 it was assumed that H, is independent of g, which is 
the case for the Green-Lagrange strain measure. This restriction, labeled (R3) in Section 
3, will be enforced below except in Section 4.5. 

4.2 Energy Variations 

As noted previously, for deriving core equations we regard the displacement gradients g 
as degrees of freedom. On substituting (4.1) and (4.3) into (4.2) we obtain the “core 
counterparts” of (3.4)-(3.6), in which v has become g: 

J=U-W = Ig r S u g + (*“-<l>) I 'g, (4.4) 

dl~C 

T= — = S r g + ^°-^ = ^-^ = 0, (4.5) 

AT = S Ag - A\£ = 0. (4.6) 

Here S , S and S denote the energy, secant and tangent core stiffness matrices, and SF 0 , 
which is independent of g, is the core counterpart of p°. 
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With this notation the first and second variations of the strain energy density can be 
expressed as 


SU = Sg T (S IJ g + <P°) + JgWg = Sg T (S r g + »°) = «g T *. (4.7) 

S 2 U = «g r S r 6g + Sg T SS r g + (6 2 gf$ = Sg T S Sg 4- (« 2 g) r *. (4.8) 

These variational equations implicitly determine S r , $ and S from S U and \&°. If the 
linearity restriction (R4) holds, the term in 6 2 g drops out as explained in the Remark 
below, and 

6 2 U = <5g r S Sg. (4.9) 


REMARK 4.1. If g = Gv with G independent of v, 6 2 g = G£ 2 v = 0 because v are inde- 
pendent variables. On the other hand, if displacement gradients are nonlinear functions of node 
displacements expressable as gi = gi{vj), then 


Sgi = = Gij Svj, 

OVj 


S'gi = 


d 2 9i 


dgi 


-SvjSvk -f — — 
dvjdvk dvj 


= Fijk Svj 6vk . 


(4.10) 


Thus Sg is still G 6v but S 2 g = (F 6v) &v, where F is a cubic array. The presence of the term S 2 g 
is taken into account in the GCCF discussed in Sections 8-10. 


4.3 Parametrized Forms 


For convenience introduce the following n g x n g matrices (with summation convention on 
i,j = 1, . . . n g implied): 

So = Eij h^hj, Si = Eij hjg Hj, Si = E{ g (h, g) Hj, ^ 

S 2 = Eij H, gg T Hji, S 2 * = En (g r H ig ) Hj, 

in which parentheses are used to emphasize the grouping of scalar quantities such as 
g r H;g. It may be then verified that, if assumptions (R3)-(R4) of Section 3 hold, the core 
stiffnesses and prestress vector in (4.6)-(4.8) possess the general form: 


S u (a,f3) = S 0 + 2 a (Si + Sf) + (1 — <a)S* + \/3S 2 + j(l — /?)S 2 4- s 0 t Hi 
= So + 2 o(Si + S^) + ( 2 — Q()Si + j/?(S 2 — S 2 ) + ^(^i 4" 

= So + ia(Si + Sf) - aSt + }0S 2 - \(1 + P)S* 2 + s&i , 

S r (<M) = So + |Si + <£Sf + (1 - <A)SJ + i(2 - tf)S 2 + }*SJ + s°i Hi 

= So + iSi + <j>Sj + (| - *)s; + }(2 - V>)S 2 + - 1)S$ + i(8°i + Si )Hi, 

= S 0 + ^Si + 4 >Sj - <f> S* + j (2 — 0)(S 2 — S*) + s,Hi, 

S = So + Si + sf + S* -(- S 2 + ^S 2 + = Sq + Si + Sj + S 2 + s,Hj, 
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Here a, (3, <p and ip are arbitrary scalar coefficients in the sense that g T S U g and S r g are 
independent of them. In fact, 

$ = S r g + *° =s { bi, (4.13) 

where b, is defined in (4.18) below. The expressions (4.12) are more general than those 
originally given by Felippa 4 because restrictions (R1)-(R2) noted in Section 3 are no longer 
enforced. Note that the secant core stiffness S r becomes symmetric if <f> — 1/2. 

The “repeatable forms” (3.1)-(3.3) of Mallet and Marcal are obtained if a = 0 = ip = 2/3 
and <p = 1/2, in which case the combinations S| + Sf + S* and S2 + |S 2 become the 
core counterparts of Nj and N2, respectively. But this observation has largely historical 
interest. More physically relevant are the following combinations: 


S d — Si + + S2, Sm = So + S/3, 

S G = St + |S 2 * + s°H, = s,H,. 


(4.14) 


These are the core versions of the initial- displacement, material and geometric stiffness, 
respectively. The core tangent stiffness is S = S 0 -f Sd + S G = Sm + S G . 

If the Generalized CCF is required for downstream element development as explained in 
Section 8, S G = s,H, is called the principal core geometric stiffness and is denoted by 
S G p. In this case the combination 


S = Sm + S G p, (4.15) 

receives the name principal core tangent stiffness. 

REMARK 4.2. Finite element practicioners may be surprised at the nonuniqueness of S u and 
S r . It appears to contradict the fact that, given two square matrices Ai and A 2 and an arbitrary 
nonzero test vector x, AiX = A 2 x for all x implies Ai = A 2 . But this is not necessarily true if 
Ai and A 2 are functions of x. More precisely, the energy core stiffness is not unique because 

g T (Si - S*)g = 0 , g T (Sf — S*)g = 0, g T (S 2 - S 2 )g = 0, (4.16) 

and the secant core stiffness is not unique because 

(Sf-sr)g = 0, (s 2 -s;)g = 0. (4.17) 

Adding u gage terms” such as those of (4.17) multiplied by arbitrary coefficients does not change 
6U and consequently the secant stiffness acquires two free parameters. Uniqueness holds for the 
tangent stiffness because the test vectors are the virtual displacement gradient variations, and S 
is not a function of 8 g. 

REMARK 4.3. Because of (4.16), an additional free parameter appears in S a if unsymmetry 
is allowed. If symmetry is enforced the first two gage expressions must be combined to read 
g T (Si +Sf -2SDg = 0. 
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4.4 Spectral Forms 

There is a more compact alternative expression of the core stiffnesses that offers theoretical 
as well as implementational advantages at the cost of some generality. Define vectors b, 
and c, as 

e« = cfg, c, = h, + |H lg , b, = ^ = hi + Hig. (418) 

Then the spectral forms (so called because of the formal similarity of equations (4.19)- 
(4.21) with the spectral decomposition of a matrix as the sum of rank-one matrices) are 


S t/ (1, 1) = s u 

= Eij c^J + s°i Hi, 

Of = /? = 1 

(4.19) 

S r (0,0) — S |^ = ^, =0 — EijhiCj +SjHi, 

(4.20) 

r, 1) = c^J + 4(s< + 4) H »i 

(4.21) 

S = Eij bibj^ + SiHi = S m + S g- 

(4.22) 


Note that S r (|,l) is symmetric but S r (0,0) is not. It is seen that for energy and secant 
stiffnesses, compactness is paid in terms of settling for specific coefficients. 


Remark 4.4. The foregoing relations may be easily verified by noting that 

Eij CicJ = So + ^ (Si + Sj ) + jS 2 , 

En bic j = So + + Sf + 2 S2 j 

Eij bibj = So + Si +Sf + S 2 , 


dg 


C i 




c. 


= St + = Eij€jHi = (Si - s°)H„ 


(4.23) 


&(cicj) _ dci_ f dc^ T _ i 


rS 2 *, 


Og 2 ” ‘ J dg \ dg 

and seeking these patterns in the general parametrized expressions (4.13). 


4.5 Generalization to H(g) 

If the H, depend on g, as it generally happens if strain measures other than Green- 
Lagrange’s are used, the secant and tangent stiffness core equations become more complex 
because of the presence of first and second g-derivatives of H^. The changes in the core 
variational equations (4.7)-(4.8) can be succintly expressed as 

8U = 6g T ((S r + S r )g + *° + ¥°) = 8g T (<& + $), (4.24) 

8 2 U = 8g T (S + S)8g + (8 2 g) T ($ + *). (4.25) 
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where S , S and <& are additional core terms that arise on account of the dependence of 
the Hi on g. 

The parametrization and efficient characterization of such terms for several strain mea- 
sures of interest in practice, notably logarithmic and midpoint strains, are presently open 
problems. Such topics would in fact be good candidates for term projects in advanced 
nonlinear finite element courses. 


5. CORE STIFFNESS EXAMPLES 


Because the core equations reflect the motion of an individual particle, their form is pri- 
marily determined by the choice of components of s, e and g that are retained in the strain 
energy density. This choice is in turn a byproduct of the mathematical idealization of the 
actual structure or structural component. 

Several cases are worked out below to illustrate the basic steps. The core expressions 
developed in these examples do not force commitment to specific elements, only to a 
mathematical model. For example the bar core equations may be subsequently used to 
develop 2 -node straight elements or 3-node curved ones. Some specific elements based on 
these equations are derived in Sections 7 , 9 and 10 . 

5.1 Bar in 3D Space 

The particle belongs to a bar moving in 3D space. The only energy contribution is due to 
the longitudinal stress. We have n<f = 3, n„ = 1 and n g = 3. To simplify node subscripting, 
Cartesian systems and displacement components will be denoted by {X,Y,Z}, {x,y,z} 
and {ux ,uy,uz] rather than {Xi,X 2 ,Xz}, {aq,:^^} and {«i,U 2 ,U 3 }, respectively. In 
the reference configuration C 0 the bar is referred to a local Cartesian system {X,Y,Z}, 
with X located along the bar axis. 

With reference to this local system, the motion of a particle initially at X is defined by 
the displacement components u x = u x (X ), uy = uy(AT) and u z = u z (X). The three 
displacement gradients that intervene in the definition of nonlinear strains are 


( 9i 


f dux/dX ' 


> = < 

duy/dX 

l 93 J 


l du z /dX J 


(5.1) 


As uniaxial strain measure we adopt the Green- Lagrange (GL) axial strain, defined as 


e = e\ = 


du 


dX 

V T 
0 
0 


+ \ 


dux 

dX 


+ 


f 9uy Y 
\dXj 


+ 


($)’ 


9 i 


9 i 


9i r + 2 { 92 


9 3 


9 3 


1 0 0 
0 1 0 
0 0 1 



— 9\ + + 92 + < 73 ) 


= h T g + yg T Hg. 


(5.2) 
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Thus for this choice of strain, = [1 0 0] and Hi = H is the 3x3 identity 

matrix. The conjugate stress measure Si = s is the second Piola-Kirchhoff (PK2) axial 
stress. The stress-strain relation is s = s° + Ee, where s° and s are PK2 axial stresses in 
the reference and current configurations, respectively, and E is Young’s modulus. 

Because H is independent of g, to form the core stiffnesses in local coordinates we can 
directly use the spectral expressions (4.19)-(4.22). First construct the vectors 


c = Cj = 


[ 1 + 2^1 ! 

r 1+01 ' 


► , b = fc>! = < 0 2 

l 703 j 

l 03 , 


(5.3) 


which inserted into the spectral forms yield 

(l + 7Sl) 2 2 ^ 2(1 +2"5l) 203(1 + !$l) 


S l/ (1,1) = £cc t + 3°H = E 


yi 


symm 


49293 

70I 


I u 

+ s 


1 0 0 
0 1 0 
0 0 1 


(5.4) 


S r (^, 1) = E cc T + s m H = E 


(l + 70 i ) 2 | 02 (l + 70 i) 703(1 + 701 ) 


2^V A » 2 * 

1 „2 
4 92 


symm 


T0203 

4*03 


+ 5 


1 0 0 
0 1 0 
0 0 1 


S = £bt> + sH = E 


(l+ 0 i ) 2 02(1 + 01 ) 03(l + 0 i) 


02 


symm 


0253 

03 



"1 

0 

O' 

+ s 

0 

1 

0 


.0 

0 

1. 


(5.5) 

(5.6) 


In equation (5.5), s m = 7(5° + s) = 3° + \Ee is the average or “half-way” stress. The 
clean separation into material and geometric (initial-stress) stiffnesses should be noted. 

5.2 Plate in Plane Stress 

As second example we consider a particle that pertains to a plate in plane stress (mem- 
brane), constrained to move in its plane. As usual we consider only the motion of the 
midplane. The Cartesian reference system and displacement components will be denoted 
by {X,Y'}, {i,y} and {«x,uy} rather than {X U X 2 }, {xi,x 2 } and {tti,u 2 }, respectively. 
The element displacement field of a generic particle originally at (X, Y) is defined by the 
two components u x - u x (X,Y) and tty = ity(X, Y). Three in-plane PK2 stresses con- 
tribute to the strain energy and four displacement gradients appear in the corresponding 
GL strain. Consequently n d = 2, n 3 = 3 and n g = 4. The four displacement gradients are 
arranged as 


g = 


/ \ 
01 


' du x /dX ' 

02 

< — J 

duy/dX 

03 

r — 1 

du x /dY 

, 04 J 


, du Y /dY , 


(5.7) 
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The strain measures chosen are the three components e; (i = 1,2,3) of the GL strains 
defined in the usual manner: 


(5.8) 









f 1 

) T 


'1 

0 

0 

o' 


ei 

= exx 

= 9i 

+ 

h(9 i 

+ 

02) = < 

0 

1 0 

1 g + ig r 

0 

0 

1 

0 

0 

0 

0 

0 

g, 








[ 0 , 



.0 

0 

0 

0. 








| 

[01 

T 


0 

0 

0 

o' 


e 2 

= e YY 

= 04 

+ 


+ 

04 ) = 

0 

0 

► g 

+ i& T 

0 

0 

0 

0 

0 

1 

0 

0 

g, 







| 

,1, 



.0 

0 

0 

1_ 










o 1 

\ T 


'0 

0 

1 

o' 

• + 

e Y x = 

92 + 

03 

+ 0103 

+ 0204 

= i 

1 

1 

f g+ig r 

0 

1 

0 

0 

0 

0 

1 

0 









.0, 



_0 

1 

0 

0_ 


(5.9) 


g, (5.10) 


from which expressions for h, and H, (i = 1,2,3) follow. For brevity, only the derivation 
of the tangent stiffness matrix will be described. Begin by forming the vectors 


bi = 


(1 + 01 ' 

( ° ) 

i 03 ) 

02 

V A 

b 2 J 0 

+ 
rH 1 

II 

cr 

S 

0 

03 

1+0! 

( 0 J 

l 1 + 04 > 

1 02 J 


(5.11) 


Then from (4.22) we get the core stiffness 

S = Eij b t bf + s.H, = S M + Sc, (5.12) 

where s, = ( i , j = 1,2,3), are the PK2 stresses in the current configuration. 

In full and using the abbreviations a\ = 1 + < 71 , a 4 = 1 + g 4 we get 


S M 


£11 a i + 2E\3d\g3 + £3353 Eh fli <72 + Ei3(d\ d^ + <72*73) + £33*2453 

E\ig\ + 2 E\ 3 d±g 2 + £ 33*24 

symm 


£12 digs + £13*21 + £2353 + £33*2153 
£l 25253 + E \3 d \ g 2 + £23*2453 + £33*21*24 
£2253 + 2£ 23 ai<73 + £33*2? 


£12*21 <14 + £l 3 * 2 l 52 + £23*2453 4- £335253 
£12*2452 + £1352 “I" £23*24 + £33*2452 
£22*2453 + £23 (*2l <24 + 5253) + £33*2152 
£22*24 + 2£ 2 3*24 52 + £3352 

(5.13) 



•Si 

0 

^3 

0 

S G = 


•Si 

0 

■S3 

0 




•S2 


symm 



S2 


(5.14) 
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5.3 Plate Bending 


This is similar to the previous example in that the structure is a flat thin plate but now 
motion in 3D space {X,Y, Z] is allowed. With this increased freedom the plate is capable 
of membrane stretching and bending. For the latter a Kirchhoff mathematical model is 
assumed. The three energy-contributing GL strains are now functions of six gradients. 
Consequently n d = 3, n s = 3 and n g = 6. The contributing gradients are arranged as 


8= \ 


The three GL strains are defined as 


e] = exx = 9l + Hdi + dl + 9z) — S 


\ 

9 1 


r dux/dX ' 

9 2 


duy/dX 

<73 
9 4 

► = < 

duz/dZ 

dux/dY 

9 5 


duy/dY 

w ^6 > 


, du z /dY ; 


rn 

0 

0 

0 

0 

10 J 


^2 = Cyy = + £5 + ffg) — \ 


(5.15) 



'1 

0 

0 

0 

0 

O' 






0 

1 

0 

0 

0 

0 





1 T 

W 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

g, 


( 5 . 16 ) 


0 

0 

0 

0 

0 

0 






.0 

0 

0 

0 

0 

0 . 






0 

0 

0 

0 

0 

O' 






0 

0 

0 

0 

0 

0 





1 T 

W 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

Si 


( 5 . 17 ) 


0 

0 

0 

0 

1 

0 






.0 

0 

0 

0 

0 

1 . 





' 0 > 

T 



'0 

0 

0 

1 

0 

O' 


1 




0 

0 

0 

0 

1 

0 


0 


g+ 

is 

T 

0 

0 

0 

0 

0 

1 


1 



1 

0 

0 

0 

0 

0 

gl 

0 




0 

1 

0 

0 

0 

0 






.0 

0 

1 

0 

0 

0 . 



£3 = e .V V + eyx = 92 + <74 + 9\9i + 9295 + 9396 — 


(5.18) 

which define h, and H t , i = 1,2, 3. When one reaches this level of bookkeeping it is more 
expedient and less error-prone to obtain the core matrices through symbolic manipulation. 
For example, the following Macsyma program forms Sjw and Sq in matrices SM and SG, 
respectively: 


hi : matrix( [1] , [0] , [0] , [0] , [0] , [0])$ 
h2: matrix( [0] , [0] , [0] , [0] , [1] , [0] )$ 
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h3 : matrix ( [0] , [1] , [0] , [1] , [0] , [0])$ 
g: matrix ( [gl] , [g2] , [g3] , [gl] , [g5] , [g6] )$ 

HH1 :matrix( [1,0, 0,0, 0,0] , [0 , 1 ,0 ,0 ,0 ,0] , [0 ,0, 1 ,0 ,0,0] , 

[0,0, 0,0, 0,0] ,[0,0, 0,0, 0,0] ,[0,0, 0,0, 0,0])$ 

HH2 : matrix ( [0 , 0 , 0 , 0 , 0 , 0] , [0,0, 0,0, 0,0] , [0,0, 0,0, 0,0] , 

[0,0,0, 1,0,0] , [0,0, 0,0, 1,0] , [0,0, 0,0,0, 1] )$ 

HH3:matrix([0, 0,0, 1,0,0] , [0,0, 0,0, 1 ,0] , [0,0, 0,0, 0,1] , 

[1,0, 0,0, 0,0] , [0,1, 0,0, 0,0] , [0,0, 1,0, 0,0])$ 
bl :hl+HHl.g$ b2:h2+HH2.g$ b3:h3+HH3.g$ 

SMrEll* bl . transpose (bl ) +E22*b2 . transpose (b2) +E33*b3 . transpose (b3) 

+ E12* (b 1 . transpose (b2) +b2 . transpose (b 1 ) ) 

+ El 3* (bl . transpose (b3)+b3 . transpose (bl) ) 

+ E23* (b2 .transpose (b3)+b3 .transpose(b2) )$ 
ratvars(g6,g5,g4 ,g3 ,g2,gl , a5 ,al ,E11 ,E12,E13,E22,E23,E33)$ 

SH: rats imp (SM)$ 

SG: rats imp (sl*HHl+s2*HH2+s3*HH3)$ 

These matrices may be automatically converted to TgK by appropriate Macsyma state- 
ments (not shown above). That output was reformatted by hand for inclusion here. For 
the core tangent stiffness this semi-automated process yields 

•Sa/( 1, 1) = £3354 + 2£i3(l + + £ll(l + <7l) 2 

Sm (1, 2) = £i 3 ((1 + 5i )(1 + 55 ) + 929a) + £3354 ( 1 + 5 s) + £ii(l4- 51)52 

5a/(1,3) = E 13 ((l + 51)56 + 5354 ) + £335456 -f £'ii(l -f 51)53 

5a/(1,4) = £2354 + £ 33(1 + 51)54 + £ 12(1 + 51)54 + £ 13(1 + 51 ) 2 

5m(1,5) = £1 2 ( ( 1 + 5i )(1 + 55 )) + £ 2354(1 + 55 ) + £335254 + £ 13(1 + 51 )52 

Sm(1,6) = £ 12(1 + 51)56 + £235456 + £335354 + £ 13(1 + 51)53 

Sm( 2,2) = £ 33(1 + 55 ) + 2£i352(l + 55 ) + £1152 

5a/(2,3) = £ 33(1 + 55)56 + £ 13(5256 + 53(1 + 5s)) + £115253 

Sm (2, 4) = £33 (( 1 + 5i)(l + 55 )) + £ 2354(1 + 55 ) + £125254 + £ 13(1 + 51 )52 

5m (2, 5) = £ 23(1 + 55 ) 2 + £ 3352(1 + 55 ) + £ 1252(1 + 5 s) + £1352 

5m (2, 6) = £ 23(1 + 5s)5e + £125256 + £ 3353(1 + 55 ) + £135253 (5.19) 

5m(3,3) = £ 335 ! + 2£i35356 + £1153 

5m(3,4) — £ 33(1 + gi)g6 + £23#4<76 + #120304 + #13 (1 + 5i )53 
5 m (3, 5) = #23(1 + 05 )06 + #330206 + #1203(1 + 05 ) + #130203 
S\f (3,6) = #2306 + #330306 + #120306 + #1303 

, 5m(4, 4) — #2204 + 2 #23 ( 1 + 01 )04 + # 33(1 + 01 ) 2 
5a/(4, 5) = #23 ( ( 1 + 01 )(1 + 05 ) + 0204 ) + # 2204(1 + 05 ) + # 33(1 + 0 i )02 
#A/(4, 6 ) = #23 ((1 + 01 )06 + 0304 ) + #22 04 06 + # 33(1 + g\ )03 
5a/(5,5) = #22(1 + 05 ) 2 + 2#2302 ( 1 + 05 ) + #3302 
Sm(5,6) = £ 22(1 + 55)56 + £ 23(5256 + 53(1 + 5 s)) + £335253 
5a/ (6, 6) = £2256 + 2£235356 + £3353 

(which can be further compacted by introducing the auxiliary symbols ai = 1 + g x and 
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a 5 = 1 + £5 as done in Section 5.3) and 


Sg = S gp = 


' 5j 0 0 

0 si 0 

0 0 

53 0 0 

0 5 3 0 

0 0 s 3 


53 0 0 ■ 

0 5 3 0 

0 0 53 

52 0 0 

0 52 0 

0 0 s 2 _ 


(5.20) 


REMARK 5.1. If the plate element to which the particle belong has (as usual) rotational free- 
doms, an additional geometric stiffness (the complementary geometric stiffness) appears in the 
transformation phase. Because of this, the core geometric stiffness (5.20) has been relabeled as 
S gp , where subscript P means “principal.” 

Remark 5.2. The core stiffness matrices may also be used for part of the formulation of thin- 
shell facet elements, with the proviso that global reference axes {X,Y,Z} are to be replaced by a 
local coordinate system {X,Y,Z} with Z normal to the element midplane. 


5.4 2D Timoshenko Beam 

Consider next an isotropic Timoshenko plane beam that moves in the (X,Y) plane. For 
notational simplicity it is assumed that the longitudinal axis of the beam is aligned with X. 
The only PK2 stresses that contribute to the strain energy are the axial stress si = sxx 
and the mean shear stress 52 = sxy- The corresponding GL strains are the axial strain 


17 


e i — e xx and the section-averaged shear strain t 2 = 7 xy = cxy + &yx ■ The constitutive 
equations are Sj = s° + Ee\ and S 2 = s° "b Ge 2, where E and G are the Young’s modulus 
and shear modulus, respectively, of the material. The treatment outlined below is slightly 
modified from that of a course term project by Alexander, de la Fuente and Haugen. 10 

The finite displacements are described in a local coordinate system that is attached to 
the initial position of the beam, as illustrated in Figure 1. Under the usual kinematic as- 
sumptions of the Timoshenko beam model (plane sections remain plane but not necessarily 
normal to the deformed centroidal axis) the coordinates of a particle in the underformed 
and deformed configurations may be written 


X = x 0 + C, 

X = x 0 + R r C x 0 


^ 0 
II 

0 

X 

O 

II 

X + Uqx ) 
UOY / ’ 

_ T COS # 

K ~ [sin# 


— sin# 
cos# ’ 


(5.21) 

(5.22) 


where uqx and uqy arc the components of the centroidal displacement vector Uo. Sub- 
tracting (5.21) from (5.22) gives the element displacement field 


u = x - X = 


f “x 1 _ f u 0 x — Y sin# 1 
1 u Y J \ ii 0 y + F(cos# — 1) J 


(5.23) 


Four displacement gradients contribute to the GL strains. Thus for this case we have 
n d = 2, n a = 2 and n g = 4. The four contributing displacement gradients are arranged in 
the usual pattern: 


{ 91 


' dux/dX ' 

92 

> — < 

duyjdX 

93 1 


dux/dY 

l 9\ J 


. du Y /dY , 


(5.24) 


For future use in Section 9 we note that the gradients can be written in terms 
section freedoms as 



9i ' 


e — Yk cos # ’ 

92 

> = i 

7 — Yk sin # 

9 3 


— sin# 

94 J 


k COS # — 1 j 


of generalized 


(5.25) 


in which e = duQx/dX is a generalized axial strain, 7 = duoY/dX a generalized shear 
strain, and k = dd/dX is the beam curvature. 

The matrix form of the GL strains is 


ei — gi + 1 + <72 ) 2 = 


r 1 ' 

T 

'1 

0 

0 

O' 

0 

> g+?g r 

0 

1 

0 

0 

0 

0 

0 

0 

0 

. 0 , 


0 

0 

0 

0 


(5.26) 
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f 01 

T 

'0 

0 

1 

O' 


1 1 

' g+|g T 

0 

0 

0 

1 


62 = 92 + <73 + <7l<73 + <7204 — S ^ 

1 

0 

0 

0 

g> 

lo. 


0 

1 

0 

0 



which define hi, h 2 , H! and H 2 . On introducing the auxiliary vectors 


' 1 +9i' 


( \ 
9 3 


r 1 + \9i ' 


H 3 1 

Q2 

0 

> , b 2 = < 

1 + 9 a 
1 + 9i 

> , Cj = < 

292 

0 

> , C 2 = < 

1 + 2 <74 1 
! + ( 

. 0 , 


i 92 J 


. 0 , 


292 ) 


the spectral core stiffness matrices and internal force vector can be written 

= Ec lC J + G c 2 cj" + s^Hi + S 2 H 2 , 

S r = EcicJ + Gc 2 c[ + + •si)Hi + ^(s 2 + s 2 )H 2 , 

S = Sa/ + Sgp, Sm = + Gb 2 b^, Sgp = sjHi + s 2 H 2 , 

$ = sibi -f s 2 b 2 . 


(5.27) 


(5.28) 


(5.29) 

(5.30) 

(5.31) 

(5.32) 


Because beam elements have rotational freedoms, a complementary geometric stiffness 
matrix appears when carrying out the transformation phase. This term is considered in 
the subsequent GCCF treatment of this element in Section 9. 

5.5 3D Timoshenko Beam: Kinematics 

The last example of derivation of core equations involve a TL 3D Timoshenko beam capable 
of arbitrarily large rotations. The following material is largely extracted from a recent 
paper by Crivelli and Felippa 7 as well as Crivelli’s thesis 5 and is continued with the DGCCF 
transformation phase in Section 10. The notation used in those references has been slightly 
edited to fit that of the present article. 

As in the 2D case, the beam is isotropically elastic with Young’s modulus E and sheax 
modulus G. The reference configuration of the beam is straight and prismatic although not 
necessarily stress free. A local reference frame n, is attached to it, with ni directed along 
the longitudinal axis (the locus of cross section centroids). Axes n 2 and n 3 are in the plane 
of the left-end cross section; these will be eventually aligned with the principal inertia axes 
to simplify some algebraic expressions. Along these axes we attach the coordinate system 
{X,Y,Z}. This description is schematically shown in Figure 2. We further define a set 
of moving frames, denoted by {a!,a 2 ,a 3 }, parametrized by the longitudinal coordinate 
X. Initially these frames coincide with {ni,n 2 ,n 3 }, and displace rigidly attached to the 
cross-sections of the moving current configuration. 

A beam particle originally at (X,Y,Z) displaces to 

x(X) = x 0 (X) + R t (X)C(Y,Z), C T = [0 Y Z], (5.33) 
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Figure 2 Kinematics of 3D Timoshenko beam element 


where x 0 describes the position of the centroid of the given cross-section, R is a 3-by-3 
orthogonal matrix function that orients the displaced cross section, and £ is a cross-section 
position vector. The displacement field is 

u = x — X = u 0 + (R T - I)C (5.34) 

where Uo(.X’) = Xo(-X) — Xo(X) is the centroidal displacement (see Figure 2). 

In the sequel 3x3 skew-symmetric matrices are consistently denoted by placing a tilde 
over their axial 3- vector symbol; for example 

0 03 —02 ( Oi I 

a= spin (a) = -o 3 0 o } , a = ^ a 2 > = axial (a). ( 5 . 35 ) 

a 2 — ai 0 J l a 3 J 

The skew-symmetric curvature matrix k is defined by k = R(dR T /dX), which is the rate 
of change of the orthogonal rotation matrix R with respect to the longitudinal coordinate. 
The curvature vector is k = axial (k). We shall also require later the variation of angular 
orientation <50, defined as the axial vector of the skew matrix R<5R r : 

<5® = R<5R T = — <5RR t , 8 © — axial (<50), (5.36) 
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All displacement gradients <7,7 appear in the GL strain measures. To maintain compactness 
the nine gradients are partitioned into three 3-vectors: 


' du x / dX 

gi = < duy/dx 
k duz/dX ) 


r du X /dY) 

g 2 = { duy/dY > , 
i duz/dY ) 


The 9-component gradient vector is g T = [gf 
directly here. Also introduce the 3-vectors 


hi = 


T 

82 


f 11 



f°l 

f°) 

0 

► , 

II 

O, h 3 = 

° 

1 0 J 



oj 1 

1 


( dux/dz ) 

g 3 = { duy/dZ \ , (5.37) 

{ du z /dZ ) 

gif], but this symbol is not used 


(5.38) 


With the help of these quantities, explicit expressions for the displacement gradient vectors 
g can be given as 

diiA 71 d Uq 

+ r t «c=-t^ + r c 


81 = 


dX ■ " dX 

T t\i_ _ roT 


(5.39) 


(5.40) 


g 2 =(R J -I)h 2 , g 3 = (R -I)h 3 . 

The only nonzero components of the GL strain tensor can be written 
ei=e n = hf gl + igfH gl , 

e 2 = 712 = 2 e J2 = hjg 1 + hfg 2 + |(gfHg 2 + 8^ h Si), 

e 3 = 7i3 = 2 e 13 = hf gl + hfg 3 + ±(gfHg 3 + gfHgj), 
where H is here the 3x3 identity matrix. Note that from the orthogonality of the rotation 
matrix R we find 

e 22 = hfg 2 + |gfg 2 

= R 22 ~ 1 + \{Rh + (R 22 - 1) 2 + R 23 ) = R 22 - 1 + |(2 - 2 R 22 ) = 0 , 

2 e 23 = hfg 3 + hfg 2 + gfg 3 

= /? 32 + R 2 3 + R 21 R 3 I + R 22 R 32 — R 32 + R 23 R 33 — R 23 = 0 ? 
and similarly e 33 = 0. This confirms that the only nonzero strains are (5.40). 

The strains (5.40) may be rewritten in a more physically suggestive form: 
ei = e n = e b + e f , 7 = 7i2+7i3, 


(5.41) 


e& = 


__ ( du 0 \ 
dX ) 


T 


h + 1 dUo 
h ‘ + 5 dX 


~T 


e 2 = 712 = 72 + 72 — ^2 ^ + ^2 C K -> (5.42) 


e f = C T K e + i/c T CC T « - C T ^e, e 3 = 713 = 73 + 73 = h 3 (f> + h 3 C* K. 


Here e^, ej are stretching and flexural normal strains, 72 and 7 3 represent bending-induced 
shear strains, and 7 2 , 73 are torsion-induced shear strains. The last term in ey represents 
a squared-curvature contribution to flexure, which can usually be neglected (cf. Remark 
9.2). The strain energy stored in the current configuration is 

U = f f Id dA dX, with IA = \Ee\ + i G (e 2 + e 3 ) + sjei + s 2 ei 2 + s 3 e 3 . (5.43) 

J L 0 J .4o 
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5.6 3D Timoshenko Beam: Core equations 

The PK 2 stresses associated with the GL strains ( 5 . 40 ) axe si = s n = s X x, s 2 = $12 = 
s X y and s 3 = 5 13 = sxz- The constitutive equations axe si = s^+Eei, s 2 = s\ + Ge 2 and 
^3 = S3 + Ge 3 . The spectral core stiffnesses can be compactly expressed in terms of the 
vectors c* = hj + ^Hgj and bj = h, + Hg^ for i = 1,2,3, where no subscript is needed in 
H = I. Applying the spectral formulas of Section 4.4 we obtain for the 9 x 9 core energy 
stiffness 


s 17 = 

rES^ 7 + G(S? + S?) 

g sy 

GS“- 


's°H 

4 h 

4 H' 


Gsf 

GSj 7 

0 

+ 

s° 2 H 

0 

0 

, (5.44) 


csr 

0 

GS"_ 


4 H 

0 

0 


where S^ 7 = c lC f, = c 2 c^, = c-jcJ’, S^ 7 = c 2 cf and S% = c 3 cf. At the residual 

level we obtain for S r a form similar to ( 5 . 44 ) except that the prestresses s° , i = 1 , 2,3 
have to be replaced by the midpoint stresses s ® +5,). The internal force vector conjugate 
to Sg is $ = S r g + 3 >° = 3*^ + $ r , in which 


<6 = < 

^ a — 1 

f aib, > 
0 

► 4> r = < 

' •s 2 b 2 + s 3 b 3 | 

1 

i 0 J 


{ s 3 b 3 J 


( 5 . 45 ) 


represent the contribution of the normal and shear stresses, respectively. 

The principal core tangent stiffness matrix S = Sm + Sgp is obtained from (4.22). The 
material stiffness is 


S m = 


E Si + G(S 2 + S 3 ) 
GSj 

GSj 


G S4 
G S x 
0 


GS5 

0 

G Si 


( 5 . 46 ) 


where S! = S 2 = b 2 b^, S 3 = b 3 bj, S 4 = b 2 b^ and S 5 = b 3 bf. The principal 

geometric stiffness is 



siH 

s 2 H 

s 3 H 


'(a? +£e x )H 

(s° + Ge 2 )H 

(s°+Ge 3 )H‘ 

Sgp = 

s 2 H 

0 

0 

= 

(•®2 + Ge 2 )H 

0 

0 


_s 3 h 

0 

0 


_(s° + Ge 3 )H 

0 

0 


( 5 . 47 ) 


T 

The contribution of (6 2 g) to the complementary geometric stiffness depends on the 
target variables in the ensuing transformation phase. Because this transformation requires 
the DGCCF, it is taken up in Section 10 after the GCCF is discussed in Section 8. 
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6. DCCF TRANSFORMATION TO PHYSICAL FREEDOMS 


The core stiffness matrices and internal-force vector given in (4.19)-(4.22) and (4.13), 
respectively, pertain to material particles of the structure. The behavior of each particle 
is expressed in terms of its displacement gradients collected in vector g. To create a 
discrete model the structure is subdivided into finite elements. Finite elements equations 
in terms of the physical DOFs collected in vector v axe constructed through a combination 
of core- to- physical transformations and integration over element domains. 

In this section we stay within the scope of the Direct CCF by assuming that the transfor- 
mations between g and v are linear. Because all subsequent developments pertain to an 
individual element, no element identifiers are used to reduce indexing clutter. 

Over an individual element the displacement field u T = (tti, u 2 , u 3 ) is interpolated as 

u = Nv, (6-1) 

where v now collects the element node-displacement degrees of freedom (DOFs) and N = 
N(X u X 2 ,X 3 ) is a matrix of shape functions independent of v. Differentiating (6.1) with 
respect to the X, and taking the first two v variations yields 

g = Gv, «5g=G6v, 6 2 g = 0, (6.2) 

(for the last one see Remark 4.1). Invariance of the strain energy variations 8U and 6 2 U 
obtained by integrating (4.7)-(4.8) over the element reference volume yields 

K u = / G T S U G dV 0 , K r = / G T S r GdV 0 , K = / G T SGdV 0 , (6.3) 

JV o J Vo JV Q 

f= / G r *dV 0 , P = f G T W dV 0 , p°= / G T *° dV 0 , (6.4) 
JVo JVo JVo 

Although the dependency of S level and ^ on g is not made implicit in these equations, 
it must be remembered that the transformation g = Gv also appears there. Because of 
the ensuing algebraic complexity, numerical integration is generally required unless the 
gradients are constant over the element. 

Often G is expressed as a chain of transformations, some of which are position dependent 
and remain inside the element integral whereas others are not and may be taken outside. 
For example, in the bar element treated below, G = T G, where G transforms g to local 
node displacements while T transforms local to global node displacements. 
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7. DCCF TRANSFORMATION EXAMPLES 


7.1 The Bar Element 


The core equations for a geometrically nonlinear TL bar were derived in Section 5.1. 
These equations are now applied to the formulation of a two-node, linear- displacement, 
prismatic TL bar element. The element has constant reference area Ao and initial length 
L 0 . The two end nodes are located at (Xi ,Y U Z X ) and (X 2 , Y 2 , Z 2 ), respectively. The node 
displacements are (t>xi» u vi> v z\ ) and (vx 2 , v Y 2 i V Z 2 )• The element displacement field in 
local coordinates { X,Y,Z } may be interpolated as 



Nx 

0 

0 


0 0 
N t 0 
0 Ni 


N 2 
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0 

N 2 
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0 

0 

n 2 



' VX1 ' 


VY 1 
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VZ 1 


V X2 
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VY 2 


, V Z2 j 


Nv, 


(7.1) 


where A r i — 1 — X / Lq and N 2 = X/L are linear shape functions. Differentiating with 
respect to the reference coordinate we get 


G 


1 

Lq 


-1 

0 

0 


0 

-1 

0 


0 

0 

-1 


1 0 
0 1 
0 0 


v=— Gv, 

-oo 


(7.2) 


This transformation may be applied to the core matrices and vectors derived in Section 
5.1. For example, application to the core tangent stiffness (5.6) yields 


K = i?S 

Jv 


G T SG dV 0 = G T (Ebb T + sH) G, 

V 0 L 0 


(7.3) 


Finally, transformation to node displacements (uxt, v Yii v Zi)> * = 1,2 is handled in the 
usual manner by writing the local-to-global transformation equation 




' 

u = < 

Uy 

f = 


l u z J 

! 


Txx T X y Txz 
Tyx T Y y Tyz 


tangent stiffness matrix in local coordinates is given by 



f UX ' 

< 

Uy 

-1 

l U Z j 

= 

1,2. 


) =Tu, 


(7.4) 





G T (Ebb T + sH)G 


T 

0 


0 

T 


(7.5) 


For this simple element all entries may be obtained in closed form and no numerical 
integration is necessary. An efficient implementation of the tangent stiffness matrix (7.5) 
in the form of a Fortran subroutine is given in Appendix 2. This implementation forms K 
with approximately 160 floating point operations. 
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7.2 Iso-P Plane Stress Element 


For the case of plane stress considered in Section 5.2, we shall asume that the associated 
finite elements are isoparametric displacement models with n nodes, and that (as usual for 
such models) the nodal freedoms are of translational type. The transformation to physical 
DOFs can then be handled within the purview of the DCCF. 

As in Section 5.1, the reference system, current system and in- plane displacement com- 
ponents are denoted by {X,Y}, {x,y} and {ux,«y}, respectively. The element nodes 
are located at }, (i = 1 , ... n) in the reference configuration Co and move to 

{xi = -Y, + u xii Vi = Yi + uyi }, ( i = 1, ... n) in the current configuration C. The 
element displacement field may be expressed as 


h* X _ 

'N , 

0 

n 2 ... 

0 ' 

l Uy J 

0 

JVi 

0 ... 

N n _ 


vxi 

vy\ 

{ V X2 > 


= Nv, 


l v Yn ) 


(7.6) 


in which iV, are appropriate isoparametric shape functions written in terms of natural 
coordinates such as f and rj for quadrilaterals. The G matrix follows upon differentiation 
with respect to X and Y, and all core equations transformed as per (6.3)-(6.4). For 
example, the physical tangent stiffness is 

K= f G t (S m + S g )G dV, (7.7) 

Jv 0 


where Sm and S a are given by (5.13) and (5.14), respectively. As in the case of linear 
elements, (7.7) is most conveniently evaluated by numerical integration. Because several of 
the integrand matrices are sparse, in the interest of efficiency in the computer implemen- 
tation the integrand may be symbolically evaluated through a computer algebra system 
such as Macsyma , Maple or Mathematica, and automatically converted to Fortran or C 
program statements before being encapsulated in the Gauss quadrature loop. 


8. THE GENERALIZED CCF 


As discussed in Section 2, the Generalized Core Congruential Formulation or GCCF is 
required when the relation between displacement gradients g and finite element degrees 
of freedom v is nonlinear. This complication occurs in elements with rotational freedoms, 
such as beams, plates and shells, if finite rotations are exactly treated. 

Recall the expression (4.8) of the second variation 8 2 U of the internal energy density. This 
expression has the core tangent stiffness S as kernel of the quadratic form in <*)g. The core 
internal force <5 also appears in the inner product (<5 2 g) 3>. This second term may either 

survive or drop out depending on the relation of g with the target physical or generalized 
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coordinates (the latter term is explained below) chosen in the CCF transformation phase. 
In the case of the DCCF, this term drops out and 

s = S m + Sc ( 8 . 1 ) 

is the tangent core stiffness, which forward transforms as per (6.3). This is the situation 
considered so fax. But if that term survives two things happen. First, (8.1) is relabeled as 

S = S\i + S gp, (8-2) 

in which S and S gp are called the principal core tangent stiffness and principal geometric 
stiffness, respectively. Second, transforming the term (<5 2 g) T $ to freedoms v produces a 
extra term in accordance with the schematics 

K = Km + Kgp + Kgc, Sm — ► Km, Sgp —*■ Kgp, (£ 2 g) T $ •-» 8v t Kgc 6v, (8.3) 


where — ♦ and t— + symbolize DCCF-transformation and GCCF-transformation-styles, re- 
spectively. As can be seen the transformation phase produces a new term Kgc called 
the complementary geometric stiffness. That term cannot be expressed in terms of the 
variation 6g of the displacement gradients. Consequently there is no “core complementary 
core geometric stiffness” S G c that can be added to (8.2). Instead it appears as a “carry 
forward term” that materializes as a quadratic-form kernel upon transforming. 

8.1 Generalized Coordinates as Generic Target 

For elements that require the GCCF treatment a one-shot transformation between g and 
v is often replaced by a multistage transformation. The degree of freedom sets used as 
intermediate targets of this process will be collectively referred to as “generalized coordi- 
nates” and identified as q. Of course the final target: element node displacements v, is a 
particular instance of such array of choices. 

In Section 2 it was noted that two variants of the GGCF, qualified as algebraic and differen- 
tial, should be distinguished in terms of consequences on the existence of physical stiffness 
equations at various variational levels. These variants are examined below. The ensuing 
development examines the transformation from displacement gradients g to a “generic 
target” set of generalized coordinates g,- collected in vector q. These coordinates are as- 
sumed to be independent , a restriction removed later. Symbols K and f are used to denote 
tangent stiffness matrices and internal force vectors, respectively, in terms of q. 

8.2 Algebraic Transformation 


The. Algebraic GCCF, or AGCCF, applies if the relation between g (source) and q (tar- 
get) is nonlinear but algebraic. We have g = g(q) or in index notation, g t = gi(qj). 
Differentiating with respect to the q, variables yields 


% = fffj’Sqj = GijSqj, or Sg = GSq, 


d 2 q da *0 (^. 4 ) 

6 9' = 'd qj dq k 6qi6qk + ~dq } = F '}kS<lMk, or S 2 g = (F <5q) <5q, 
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Here (F6q) is the matrix Fijk 8qk = Fkij 8qk\ F being a cubic array. The array G receives 
the name tangent transformation matrix. The second term in the expansion of 8 2 gi vanishes 
because the q , are assumed to be independent target variables. 

Enforcing invariance of 8 2 U yields the tangent stiffness transformation 


K = J |g t (Sm + S gp)G + Q j dVo = Km + K gp + Kcc = Km + Kg, 


(8.5) 


where the entries of Q are (cf. Remark 4.1) Qij = Qji = Fkij^k with summation on 
k = rig. Note that Q is symmetric because Fkij = Ffcy, . Integration of Q over Vo 

yields the complementary portion Kgc of the geometric stiffness Kg. 

The internal, applied and prestress force vectors transform according to the formulas in 
(6.4) with the G defined in (8.4): 


f= f G t $ dV 0 , p= [ G T V dV 0 , P°= / G T *° dV Q . 
JV o JVo JVo 


(8.6) 


What happens to and K r ? They can be obtained, somewhat artificially, by construct- 
ing the matrix equation 

g = Wq, (8.7) 

where W is called a secant transformation matrix. Generally this matrix is far from unique 
because its n g x n q entries must satisfy only n g conditions. (Care has often to be given to 
the qj — ► 0 if 0/0 limits appear in W.) Using (8.7) we can proceed to form 


K 


u 


= f W T S u WdVo, 

Jv o 


K r = f ( 

JV , o 


G r S r W dV 0 . 


( 8 . 8 ) 


Because in general W ^ G, symmetry in the secant stiffness K r cannot be expected even 
if S r is symmetric. 

REMARK 8.1. The AGGCF is applicable to finite elements with degrees of freedoms that in- 
clude fixed- axis rotations , because such rotations are integrable. Examples are provided by two- 
dimensional beams as well as plane stress (membrane) elements with drilling freedoms if only 
in-plane motions are allowed. 

REMARK 8.2. Why is Kgc called a geometric stiffness? Because it vanishes if the current 
configuration is stress free, in which case the core internal force $ vanishes and so does Q. 

8.3 Differential Transformation 

The Differential GCCF, or DGCCF, is required if the relation between g (source) and q 
(target) is only available as a non-integrable differential form between their variations: 


8g, = G l} 8q } , or 8g = G (5q, 


8 2 g, = 


dGjj 
dq k 


8q j 8q k = F ijk 8qj 8q k , or 8 2 g = (F<5q) 6q. 


(8.9) 
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The transformation equation (8.5) still applies for K whereas (8.6) holds for the force vec- 
tors. But no integral g = g(q) as in the AGGCF exists. Consequently K u and K r , which 
require a secant matrix relation of the form (8.7), cannot be constructed. Furthermore 
Q is not necessarily symmetric; a condition for that being Fkij = F^ji or equivalently 
dGki/dqj — dGkj/dqi. 

REMARK 8.3. For mechanical finite elements the DGCCF naturally arises when three- 
dimensional finite rotations are present as nodal degrees of freedom, because such rotations are 
non-integrable. 

REMARK 8.4. The relations (8.9) have points of resemblance with the case of non-holonomic 
constraints in analytical dynamics. 

8.4 Multistage Transformation 

Up to this point the q have been assumed to be independent variables. But as previously 
noted, for complicated elements the GCCF transformations axe more conveniently applied 
in stages. The target variables in one stage become the source variables for the next one. 

What happens if the q are intermediate variables in a transformation chain? If the q 
are linear in the final independent degrees of freedom v, all previous formulas hold be- 
cause the DCCF applies for the remaining transformations, which are strictly congruential. 
But if the q are nonlinear in v, or only a non-integrable differential relation exists, term 
( dgi/dqj)S 2 qj = Gij 8 2 q } in the second of (8.4) survives. The net effect is that the geo- 
metric stiffness acquires a higher order component, implicitly defined as the kernel of 

$ t G ij 8 2 q j dV 0 , (8.10) 

This term cannot be resolved ( “resolution” meaning explicit extraction of its stiffness kernel 
in the form of a complementary geometric stiffness) until the transformation chain reaches 
downstream variables that either are the final degrees of freedom (and thus independent), 
or depend linearly on such. It is difficult to state detailed rules that encompass all possible 
situations. Instead the treatment of the 2D and 3D beam element transformations in 
Sections 9-10 illustrates the basic techniques for “carrying forward” terms such as (8.10). 

9. A 2-NODE 2D TIMOSHENKO BEAM ELEMENT 


We continue here with the derivation of a 2D, isotropic Timoshenko beam element started 
in Section 5.4. This example serves to illustrate the Algebraic GCCF. The specific element 
constructed here has two end nodes, six degrees of freedom, and reference length Lq. The 
cross section area A = A 0 and moment of inertia I = J A Y 2 dA are constant along the 
element. Axis X is made to pass through the centroid so that f A Y dA = 0. Furthermore 
it is assumed that the cross section is doubly symmetric so that f Y 3 dA = 0. 
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The element displacement field, defined by «o y(X) and #(X), is interpolated with 

linear shape functions: 
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where N\ = 1 — (X/Lq) and iV 2 = 1 — N x = X/Lq. Consequently 
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9u 0 x 


V X2 - 


ax T 0 

are constant over the element. 
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VX2 ~ VX± 

Lo 
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dO O 2 — 0\ 


dX 


To 


(9.1) 


(9.2) 


9.1 Generalized Coordinates and Stress Resultants 

As intermediate set of generalized coordinates we take q r = [e 7 k 0\. These four 
quantities are constant over each cross section and may be viewed as cross-section orien- 
tation coordinates. Consequently when obtaining stiffness matrices and internal forces in 
terms of q it is convenient to integrate over the beam cross section. The resulting quanti- 
ties appear naturally in terms of cross section stress-resultants as shown below. In terms 
of these generalized coordinates the auxiliary vectors bj listed in (5.28) become 


b, 


( 1 + 91 ' 


1 + e — Y k cos 9 y 

92 

> ~ < 

7 — Yk sin 0 
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0 ’ 

l 0 j 
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9 3 

1 + <74 
1 + 01 

92 


-sin 6 
cos 0 

1 + e — Yk cos 0 ' ’ 

7 — Yk sin# ; 

(9.3) 


The well known stress resultants of beam theory are the axial force N, transverse shear 
force V and bending moment M. They are obtained by integrating the PK 2 stresses over 
the beam cross section: 


N = 
V = 
M = 



s l dA = EA(e + i(e 2 + 7 2 )) + \EIk 2 + N° , 
s 2 dA = GA a u^ + V°, 
siY dA = — EIku> ( + M°, 


(9.4) 


where u ( = (1 + e)cos# + 7 sin# and = 7 C 0 S# — (1 + e)sin# can be viewed as 
generalized skew strains. In (9.4) IV 0 , V° and M° denote initial-stress resultants (stress 
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resultants in Co, also called prestress forces), A = A 0 , I = J A Y 2 dA, and A s = // A , in 
which /x is the usual shear correction factor of Timoshenko beam theory. Because of the 
doubly-symmetric cross-section assumption, a term containing the third-section-moment 
J A y 3 dA has been omitted from the expression for M. 

In addition to N , V and M, the following higher order moment, which is absent from the 
linear theory, appears in the residual force and tangent stiffness: 

C = [ s x Y 2 dA = El((e + I(e 2 + 7 2 )) + \EHk 2 ) + C°, (9.5) 

J Aq 

in which 7i = f A Y 4 dA. If terms in k 2 are neglected, 

C - C° = (N - N°)(I/A) = (N - N°)r 2 , (9.6) 

where r = yJIfA is the radius of gyration of the cross section. If such terms are retained 
this relation is only exact if r 2 ='HjI and approximate otherwise. 

Remark 9.1. One may verify that f s^Y dA vanishes identically. This serves as a check of the 
strain distribution equations. 


9.2 Transformation Matrices 


The differential relations required to establish the tangent transformation are obtained 
from (5.25) as 
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Gi <5q, 


G 2 6v, 


The transformation relating 8g = G 8v may be obtained as the product 


(9.7) 


(9.8) 


G = GiG 2 



0 

-1 

0 

0 


T(cos 9 + (L 0 — X)k sin 9) 1 0 

F(sin 6 — (Lq — X )k cos 0) 0 1 

—(Lo — X)cos8 0 0 

— (Lq — X)sin# 0 0 


Y( — cos 9 + X k sin 6 ) 
Y( — sin ^ — Xk cos 0) 
—X cos 9 
—X sin 8 

(9.9) 


but it is more instructive (as well as conducive to higher efficiency in the computer imple- 
mentation) to perform the transformation phase in two stages. 

Observe that the first transformation (from g to q) is nonlinear and algebraic whereas the 
second one (from q to v) is linear. Consequently we have to use the AGCCF for the first 
transformation but the second one can be done simply through the DCCF. 


30 



9.3 Internal Force Vector 


The internal force vector in terms of q, denoted by f 9 , is obtained from the core expression 
(5.32) for and the matrix Gi given in (9.7): 


f 9 


= JgT 

Ja o 


[ ft] 


' N( 1 + e) — M/ccos# — V sin# 1 
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Finally, application of (9.8) and integration over the element length yields 
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This vector satisfies translational equilibrium. 


(9.10) 


(9.11) 


9.4 Tangent Stiffness Matrix 

Transforming to generalized coordinates q produces three components of the tangent stiff- 
ness matrix: 

K 9 = f (G^(Sm 4- Sg)Gi + Q) dA = + K^p + (9-12) 

h i 0 

The entries of K 9 W , obtained through symbolic manipulation, are 


1) = EA{ 1 + c ) 2 + GA, sin 2 9 + EIk 2 cos 2 9, 

K q M { 1,2) = EA{ 1 + t)~i — GA, sin 9 cos 9 + EIk 2 sin 9 cos 9, 

I< q M {l,3) = EIk ((1 + e)(l + cos 2 9) + 7sin 9 cos 9), 

K\j( 1,4) = EIk 2 u-, cos 9 + GA,u> e sin#, 

Ii q M { 2,2) = EA~f 2 + GA, cos 2 9 + Eltsin 9, 

K q M (2,3) = EIk ((1 + e) sin 9 cos 9 + 7(1 + sin 2 0)), 

Km( 2, 4) = EIk 2 u-y sin 9 — GA s u t cos 9, 

Ii q M ( 3,3) = EIuj 2 + EHk, 

K q M {3,4) = EIk ((1 + c) 7 (cos 2 9 — sin 2 9) + (t 2 - ( 1 + e) 2 ) sin 9 cos 9 ) , 
K q M (4,4) = EIk 2 4> 2 9 + GA 8 u ( 2 . 

The principal geometric stiffness, which is readily worked out by hand, is 
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(9.14) 
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K 9 — 
"• CP ~ 
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symm 


— M cos 9 M k sin 9 — V cos 6 
—M sinO —MncosO — VsinO 
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The new term contributed by the AGCCF to K 9 is the complementary geometric stiffness 
K g GC . Its source is the matrix Q introduced in Section 8.2. The entries of Q are Q tj = 
( d 2 9 k/dqidqj)$ k , where the components of g and # = s 1 b 1 + s 2 b 2 may be obtained from 
(5.25) and (9.3), respectively. 

The entries of Q were symbolically generated by the following Mathematica module: 


Qmatrix0f2DTimoBeamElement [eps_ ,gamma_ ,kappa_ , theta. ,Em_ ,Gm_ ,Y_] : = 

Module [-[g ,hl ,h2 ,H1 ,H2,el ,e2,sl ,s2 ,bl ,b2,phi , i ,j ,k>, 
q= -Ceps .gamma, kappa, thet ah phi=-Cl , 1 , 1 , 1>; 
g={eps-Y*kappa*Cos [theta] , gamma- Y*kappa€ in [theta] , 

-Sinftheta] , Cos [theta] -1} ; 
gg=«g[[l]]>.{g[[2]]>,{g[|3]]>,{g[[4]]»; 
hl={{l> , {0> , {0} , {0» ; h2=«0> , {1> , {1> , -[0» ; 
H1=-C-C1,0,0,0>,{0,1,0,0>,{0,0,0,0>,{0,0,0,0»; 
H2={{0,0,1,0>,{0,0,0,1>,-[1,0,0,0>,{0,1,0,0»; 
el=(Transpose [hi] . gg+ ( 1/2) *Transpose [gg] .Hl.gg) [[1,1]] ; 
e2= (Transpose [h2] .gg+(l/2)*Transpose[gg] .H2.gg) [[1,1]] ; 
sl=Simplify [Em*el] ; s2=Simplify [Gm*e2] ; 

bl={l+eps-Y*kappa*Cos [theta] , gamma- Y*kappa*Sin [theta] ,0,0}; 
b2={-Sin [theta] ,Cos [theta] , l+eps-Y*kappa*Cos [theta] , 
gamma-Y*kappa*Sin [theta] } ; 
phi=Simplify [sl*bl+s2*b2] ; 
q=Table[0,{4},{4}]; 

For [i=l , i<=4 , i++ , For[j*l, j<=4, j++, For[k=l,k<=4,k++, 

Q[[i,j]]»Q[[i,j]] + (D[D[g[[k]],q[[i]]],q[[j]]])*phi[[k]]]]]; 

Return [Q] 

]; 

The output of this module was integrated over the cross section and pattern matched with 
the expression of the stress resultants (9.4)-(9.5) to produce 


0 


Khc = 


symm 


0 0 0 

0 0 0 

0 —Mu-f 

— Vu> 1 + M KU) f 




Cl c 2 


(9.15) 


which added to (9.14) yields the geometric stiffness 


K 


7 __ 
G — 


N 0 
N 


symm 


—M cos 9 Mk sin 0 — V cos 9 
— M sin 9 — Mk cos 6 — V sin 9 

C -Mu> 7 

V UJy -f~ Af KU) € 


(9.16) 
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Finally, the tangent stiffness in terms of q is K ? = K m + Kg- Denoting the entries of K ? 
by Klj, i,j = 1, ... 4 the tangent stiffness matrix K in terms of node displacements v is 
formed through the DCCF transformation 


K 


-l 


L o 


GTK q G 2 dX = 


I\ g 

IX ll 


l<n 

A& 


A?, - lL 0 K’ t 

Kh - \UKl, 

K! , - L,K !, + ±LlKl 


symm 


~K q n 

-K’ n 

-A?, - \U K\, 

-K\2 

-A? 2 

-K’ 3 - \UKl, 

, + 

-iq 3 + i L„KI 

— A?, + iLlKt, 

K'u 

I< q 

IX 12 

K< 3 + 1L„A'? 4 


K q 

1X 22 

A?, + 1IoA? 4 


k;, + L„iq t + J 


(9.17) 


The above rule can be applied to and Kg should separate formation be desirable, as 
when setting up a stability eigenproblem. Using these schemes K can be formed at a cost 
of approximately 300 floating-point operations per element, which is not too different from 
the cost of a TL 3D bar. 


If the reference configuration is not aligned with X, the preceding expressions apply to the 
local system {X, Y }. A final local- to-global transformation step, similar to that discussed 
for the 3D bar in Section 7.1, is then necessary. This step can be handled by a simple DCCF 
transformation, because the finite rotation 0 remains the same in global coordinates. 


REMARK 9.2. The foregoing exact expressions contain curvature-squared terms typically in the 
combination Ik 2 . This can be shown to be of order ( r/R ) 2 compared to other terms, where r 
is the radius of gyration of the cross section and R = 1/k the radius of curvature of the current 
configuration. For typical beams {r/R) 2 is 10“ 6 or less; consequently all such tiny terms may be 
dropped without visible loss of accuracy. For highly-bent extremely-thin beams, however, that 
ratio may go up to 0.01 in which case the k 2 terms might have a noticeable though small effect if 
retained. 


9.5 Can a Secant Stiffness be Constructed? 


To attempt the construction of a secant stiffness K rq in terms of generalized coordinates q 
one should obtain a secant matrix form of the relationship g = g(q). As noted in Section 
8.2, such form is far from unique. One possible choice is 


{ 91 ) 

'1 

0 

—Y cos 6 

0 

\ 92 1 _ 

0 

1 

—Y sin 9 

0 

u r 

0 

0 

0 

— sin 6/9 

l 94 . J 

0 

0 

0 

(cos B — l)/0 _ 


6 


7 

K 

6 


w iq , 


(9.18) 
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Table 1. Internal energy and its variations for 3D Timoshenko beam element 


Core 

Particle 

g 

Section Gradients 

Cross-Section 

w 

Section Orientation 

Cross-Section 

z 

Physical DOF 
Whole Element 

V 

U = ig T S l/ g + g T *° 

— 

— 

— 

6U = 6 g T (S r g + ¥°) 

£ 

II 

6U Z = 6z T f z 

SU = Sv T f 

I 

8 2 U = 6g T S <5g + 6 2 g T # 

6 2 Ug = Syt T S 6w + T 

S 2 U Z = 6z t K z Sz 

6 2 U = 6v t K6v 


which has the merit of not being too dissimilar from Gj. Note that some care must be taken 
as regards some 0/0 limits. Then K r? = J A G^S r Wj dA, which may be easily worked out 
in closed form but is unsymmetric. Because q is linear in v, the next transformation is 
simply K r — f 0 ° G^K r? G 2 dX which can be handled through a scheme similar to (9.17) 
but with an unsymmetric kernel matrix. 

10. A 2-NODE 3D TIMOSHENKO BEAM ELEMENT 

We continue here the development of a two-node 3D Timoshenko beam element started in 
Section 5.5. As can be surmised, the development is more complex and demanding than 
for its 2D counterpart. Only a summary taken from Crivelli’s thesis 5 and Crivelli and 
Felippa' is presented here. The transformation phase to pass from the core equations to 
the element nodal degrees of freedom is carried out in three stages: 

1. From particle displacement gradients g to generalized gradients w at each cross sec- 
tion. An integration over the cross section area is involved. 

2. From generalized gradients w to cross-section orientation coordinates q. The rota- 
tional parametrization is introduced at this stage. 

3. From cross-section orientation to finite-element nodal degrees of freedom v. An inte- 
gration over the element length, as defined by the shape functions, is involved. 

These transformation stages are summarized in Tables 1 and 2, which together also serve 
to define notation 

10.1 Transformation to Generalized Gradients 

The first set of target variables are the generalized gradients w(X) at each reference cross 
section defined by the longitudinal coordinate X. The components of w are indirectly 
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Table 2. Core-to-physical-DOFs transformations for 3D beam element 


Core Level 
Particle 
g 

Section Gradients 

Cross-Section 

w 

Section Orientation 

Cross-Section 

z 

Physical DOF 
Whole Element 

V 

$ 

n= j w T *dA 

f, = z T n 

f L ° r 

f= / Gji z dX 

S 

J A o 

S= I W T SW dA 

K : = Zi T S% + SgCz 

Jo 

K = [ ° GjK.G, dX 


J Aq ! 

' 

1 

Jo 


Sg = W 8w 6vs = Z 6z Sz = G 2 <5v 


( 10 . 1 ) 


given through their first variation: 

6w = {^ S1 IDT * 0 } ’ (10,1) 

where 80, defined in (5.36), measures the variation of angular orientation. Because this 
quantity is not generally integrable for three-dimensional motions, it is not possible to 
express © as a unique function of the displacements. The variation of g 3 is 


eg, = ^ + r t c t ^ + R t C T *e& + R T Ji0c T K, 

dX aX 


( 10 . 2 ) 


where we used the relation 5 8k = d80/dX + k80. On using the commutative law ab = 
b a and Jacobi’s identity ab = ab - ba we may rewrite (10.2) as 


. d d)Un t’-T 9 8® ■nT~T'j-c^~\ 

igl = ~§F + R c !F + c 


(10.3) 


For the other gradient vectors we have Sg 2 = ^R 7 ^ = R T ^0h2 = R h 2 <5© and Sg 3 = 


~ T 

R T h 3 80, which can be collected in matrix form as 


r <tei 

= l 


1 r t c t r t « t c ( Tot I Wl 

0 0 R T h^ < <±80 > = w 2 8w = W <5w, (10.4) 


0 0 


rr~T 

R T h, 


d<5uo 

~1T 


where I is the 3-by-3 identity matrix and Wj are 3-by-9 matrices. The second variation of 
g, which is required for the complementary geometric stiffness, is 




S 2 g, = R' 6©< 


180 rp- T ~T dSO 

i^ + R T d©C k80 + R T 80C -TT7- 
dX dX 


+ R T 80~C k80 + 8 2 R r k + R t <? 8 2 k, 


(10.5) 


8' 2 g 2 =8 2 R T i 2 , 8 2 g 3 — 8 2 R T i 3 
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At this point it is appropriate to introduce the following section resultants: 


V=Aa b + T°, 

«S b E } 


Q = fi s A+Q°, 

t=t 2 + t 3 , r 2 = G7 2 h 2 , r 3 = G7 3 h 3 , 


M,=EI s K e +Ml, 

I5 = / CC dA, K e — cj) 

Ja 0 

(10.6) 

Ad T = fi t GIpK + Ad”, 

lr = f ifiA. 

Ja« 



Here V , Q, Ad^ and Ad r are axial forces, shear forces, bending moments and torsional 
moments, respectively, at the current configuration C; "P”, Q°, Ad” and Ad” axe similar 
quantities at the reference configuration Co', fi a and fit are transverse- shear and torsion 
coefficients that account for the actual shear stress distributions, respectively; and Is and 
Ip are the cartesian and polar inertia tensors, respectively, of the cross section. Should 
the axes Y and Z be aligned with the principal inertia axes the latter simplified to 



'0 

0 

0 ' 


CO 

+ 

<N 

0 

0 ' 

I s = 

0 

I 22 

0 

> IP = 

0 

hi 

0 


0 

0 

1 

co 


0 

0 

I 22 _ 


(10.7) 


Because the relation between g and w is of differential type the applicable transformation 
rules are those the DGCCF, and no energy or secant stiffness survives. Thus only the 
internal force vector H and tangent stiffness S associated with w are derived below. 

Internal Force Vector. The generalized internal force vector is 

n= f w T & dA = V] [ si w T b, dA = n a + n T , ( 10 . 8 ) 

Ja q t J Aq 


where lZ a and TZ r are the contributions of the normal and shear stresses respectively. 
Detailed calculations result 5 in the following exact expressions: 



-R T (V<f> + kM a y 


R t Q 

n„ = 

4> T M a 

K T e M a 

, Fi T = 

M r 

_4> Q + k T Mr _ 


(10.9) 


For small deformations in which the squared curvature may be neglected, R « I, <j> hi, 
K e ~ k and itA d^ ss 0. If these approximations are made, 



~T 


Q 

n, = 

h, M a 
0 

II 

k 

£ 

Ad r 

~T 


( 10 . 10 ) 
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These resemble the classic linearized theory equations. Furthermore observe that the term 
VR T <t> corresponds to the internal force of the TL 3D bar. 

Tangerit Stiffness. For the tangent stiffness we have the decomposition 

S = Sm + Sgp + &GC- (10.11) 

Furthermore, since w is nonlinear in downstream variables, the complementary geometric 
stiffness splits into two components: 


SgC = $GCw + $GCq, (10.12) 

where Sgcw and Sccq contains terms that depend on the first and second variations, 
respectively, of R and k. The notation is suggested by the fact that Sgcw can be merged 
into Sgp to yield the geometric stiffness Sg w = Sgp + Sgcw, which is associated with 
the generalized gradients w and independent of the rotational parametrization selected in 
the next set of target variables q. On the other hand, the kernel S G Cq cannot be extracted 
at the w level and must be carried forward to the q level because it is parametrization 
dependent. Each of the components in ( 10. 1 1 )-( 10. 12) may be expressed as the sum of two 
contributions, one from the normal stresses and one from the shear stresses: 


«Sjv/ — Sm<j + Sm-t, Sgp — Sgp^ + Scp r , Sgcz = Sgcxo + Sgcxt , x = w,q. 

(10.13) 

Material Stiffness. The generalized core material stiffness is given by the congruential 
transformation 


Sm = 



W r S M W dA 


J2 f E l W T b l bJWdA = S M ,+S Mr . 

i d Ao 


(10.14) 


Carrying out the algebraic manipulations one obtains 


\R T (4><f> T + k T I s k)R 


Smct = E 


R T kl s 4> 

- T 

<f> I S(f> 


symm 


R T klsK-e 

~T ~ 

^ Is^Ce ? 

~T t ~ 

K e I S K e J 


(10.15) 



'AR t I_lR 0 

AR T Ij_<p 


'0 

0 

O' 

Smt — t*G 

Ip 

Ipk 

, in which lx = 

0 

1 

0 


symm A<j> I±<j) + K T I P K 


0 

0 

1_ 


(iO- 16 ) 

The contribution R T 0<£ T R is the core material stiffness of a TL 3D bar. 


Geometric Stiffness due to Normal Stresses. It is convenient to work out together all 
geometric stiffness terms produced by the normal stresses, i.e. 


Sg<j — Sgp, t + ScCwa + ScCqa = Sgw, 7 + ScCqcr- 


(10.17) 
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The appropriate definitions are 


>GP<r 


-L 


suWfHWj^, 


$gc<t= / ^nbi 8 2 g dA = 8w T S GCw<r 8w + V(8 2 R,8 2 k), 

J Ao 


(10.18) 


where T contains S GCq as q level kernel. Carrying out the algebraic manipulations one 
obtains 

T 


^Gw(t = $GP<t + SgCuht = 


VI R/ M. 


R T k T M, 


(10.19) 


0 _ _ 

L symm <j>M<,k + k T J 

The term VI corresponds to the core geometric stiffness of the 3D T1 bar. 

The higher order term in (10.18) may be expressed as 

V a {8 2 R,8 2 K ) = M^4 >S 2 k + <t> T R8 2 R T kMvSq 1 ' (v(j> T M a ) + U(kM»; <£)) 8q, 

_ ( 10 . 20 ) 

Consequently 

Sccqv = V(4> T M a ) + V{kM (r \ <f>). (10.21) 

Because the next-level target variables q include the finite rotation parametrization, matri- 
ces V and U depend on that choice. They are the source of unsymmetries in the stiffness 
matrices when certain rotational parametrizations axe adopted, such as the incremental 
rotation vector. If the rotational vector is chosen these matrices are symmetric. 

Geometric Stiffness due to Shear Stresses. The contribution of the shear stresses to the 
geometric stiffness is 


^ Gt &GPr “H ^GCwt + &GCqr — &Gwr GCqr • 

The appropriate definitions are 

Sgpt= / 5 12 (WfHW 2 + W 2 HWi) + s 1 3(WfHW3 + W 3 HW 1 )d^ 

J v4o 

Sgct= / (si 2 b 2 +5 13 b 3 )6 2 g dA = 8w t S GCwt Sw + T t (8 2 R,6 2 k). 

J Ao 


( 10 . 22 ) 


(10.23) 


Carrying out manipulations one obtains the surprisingly simple form for S Gt 

[ 0 0 R t Q T 

$Gwt = S G p r + SqCwt = 0 0 

. symm 0 

The terms due to the second variation of g become 

Tr = Q T 8 2 R$ + M?8 2 k. 

The kernel carried forward to the q level is 

«S GC?r = V(jVl r ) + U(Q;*). 


(10.24) 

(10.25) 

(10.26) 
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10.2 Transformation to the Rotational Vector 

The second transformation stage passes from w to z, which is a vector of generalized 
displacements, also associated with a beam section, which embodies the parametnzation 
of the cross section rotation: 


z — 


duo da 




r f dS\io dSa \ 

Sz = \ 1 x lx 4 “J ■ 


(10.27) 


dX dX 

Here a denotes the rotational vector parametrization defined by the standard formulas 

a = axial (a), R = exp(d T ), (10.28) 

and which may be extracted from R by 

arCSin(r) a*ia l (R T -R), r = i||axial(R r - R)||. (10.29) 


ac - log R = 


2r 


Because only the variations of w are known the relation between w and z is also of 
differential type: 


Sw = ZSz, or <5w = 


10 0 

o Y ( z) qp 

0 0 Y(z) J 


1 

d8 Uo i 

< 

d8a 

■ajr 


l Set J 


(10.30) 


in which 


sin a 


in|a|\ao! T 1 — cos |a| , 


or 


( sinlo!l\ . . 

Y(a) = -r^i + (i--^) w — 


a 


(10.31) 


On applying the transformations (10.30) we find for the internal force and the material 
and principal-geometric components of the tangent stiffness matrix. 


f = Z T (Ha + Kt), Km, = Z T («S M )Z, K G p, = Z T (5 Gu) Z. 


(10.32) 


The materialization of the geometric stiffness terms Sccqa an d ^GCgr for the rotational 
vector needs additional work. We state here only the final result: 


U(r;$) = 


1 

O 

o 

0 

1 


l 

o 

o 

o 

0 0 

T(M r ) = 

0 T[ 

symm U T 


. symm T£ . 


(10.33) 
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where 


U r = + c 2 + $> T rj 4- c 3 + r^d^I + <&a T + a4> T T^ 

4- c 5 (a T $> ( ra T 4 - ar T ) + a T r^a$ T 4- <&a T ^ + r^aa 7 ^^ 

4- c 4 /cii$aa T + cqt^ aa T $> aa T , 

V[ = c 2 A"l r 4- csaAd^ 4- c 5 aa T A4 r 4- c 7 (A4 r a T 4- a T Al r l) 4- c 8 a r A / t 7 -aa T , 


da 7 da r 

V5 = _ C3 _ M.I-c 4 - 


t I t-» dot y y dot 

X r aa 7 +c 5 [ — M r a T + a— M.+a 7 — A 4 r I 


\ 

T’da T /da . , r . . da T \ 

+ c,a ^M r ca + c,^—M T +M rlI j + 

+ C, 3? o(aM' + M r o I ' + o , ’Mj) + 


rp 

dot T Kjt T 

c«— aa M r aa , 

dX 


in which 


(10.34) 


Cl = 

sin a 

C2 = 

1 — cos a 

sin a — a cos a 

* 

a 

c \ + 3 c 3 

7 ’ 

cr 

Ci + 2 c 2 

£3 0 i 

C3 + 4C5 

c 4 = 

9 > 

a z 

1 + C\ 

c 5 = 

7 ’ 

a z 

3c 3 — 2c 2 

^6 2 ’ 
a z 

C 5 - 5 c 8 

C7 = 

9 1 

a 1 

c 8 = 

2 ’ 
a 1 

c 9 = 2 • 

a* 


(10.35) 


A similar approach can be taken with (10.30), which defines T a . The tangent stiffness 
matrix can be obtained by superposing all contributions. 

10.3 Transformation to Finite Element Freedoms 


The final stage introduces a finite element representation for the degrees of freedom. The 
beam or beam assembly is divided into a set of two-node finite elements. Each of these 
nodes has three displacement degrees of freedom and three rotational degrees of freedom 
corresponding to the three {ax,ay,ftz} components of the rotational vector a. Each 
element in turn has twelve freedoms which are collected in the array \ T = {u n ot n } T 
where d n collects the six translational freedoms while a„ collects the six rotations. The 
cross-section state vector z is approximated inside each element by 


z = 


N 0 

n <*N 

0 JX 

0 N 




v. 


(10.36) 


where N is a matrix of linear shape functions. Since 6q = G z <5v the final internal force 
vector f and tangent stiffness matrix K of each element are obtained through the DCCF 


40 



transformations 


f = 




G^K z G z dX. 


(10.37) 


The choice of shape functions for the rotational vector poses some subtle questions. In 
small-deflection analysis it is common practice to select all Timoshenko beam shape func- 
tions to be linear in X. This choice obviously enforces nodal compatibility while pre- 
serving constant curvature states. But for finite deflections a linear interpolation for the 
rotational vector components cannot exactly represent a constant curvature state unless 
the rotations are about a single axis (plane rotations). The same is true if the rotation 
matrix RfX) is interpolated linearly. On the other hand, linear interpolation of Euler 
parameters does preserve the constant curvature state. This motivated the development 
of an interpolation scheme that starts from the 4 Euler parameters Ci(-Y), i — 0,1,2, 3, 
€ i = 1 that orient the normal of a cross section at X. These are collected in the 4- vector 
e = { e 0 €3 } . Given the eight end values e(0) and e(L) the interpolation that 

can copy a constant curvature vector k is found to be 5 


e(C) = cos(C) 



tan(£) 

tan(Ci) 



+ 


shi(C) 

sin(Ci) 


<L), 


(10.38) 


where C = ^kX , k = v k t k. The constant curvature vector can be extracted 

from the end values through the formula 


* = jpi [( e W - 2eo(£)l) c(0) - (e(0) - 2e o (0)l) 6(0)] , (10.39) 

This interpolation is then transformed to the variations in terms of the rotational vector. 
Details are provided in Reference 5 . 


11. APPLICATION EXAMPLES 


Several application examples solved with the Timoshenko beam elements are described 
below to show that they do not suffer from the restriction to moderate rotations that 
several authors attribute the TL description. 

11.1 Cantilevered Beam under End Moment 

This is a classic validation test for geometrically nonlinear beam and shell elements. A 
cantilevered beam of initial length Lq is forced into pure bending by application of an end 
moment M. The beam bends into an arc of circle with curvature k = M/(EI) and end 
rotation 9 en d — MLq/{EI). The test were run with E = 30.0 10 6 , G = Ej 2 , A = 1 , 
I = 1/12, Lq = 1000, and M = 15708. 4A, where A is a load parameter. The {X, Y] axes 
are placed at the left end of the beam. The load scaling is chosen so that for A = 1 the beam 
bends to a full circle of radius Lq/(2tt) = 318.31. The results for the tip deflection obtained 
using 10 elements along the length are shown in Table 3. The results for the 2 D elements 10 
and 3 D elements 5 were essentially the same within solution acceptance tolerances. The 
average number of full-Newton iterations per step was 4 . 
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Table 3 Computed solutions for plane cantilever beam under 


pure moment 


Load level 

Numerical solution 

Analytical solution 

A 

X Up 

Y tip 

®tip 

x itp 

^ tip @tip 

0.25 

-362.74 

637.29 

1.57 

-363.38 

636.62 1.57 

0.50 

-1000.03 

639.23 

3.14 

-1000.00 

636.62 3.14 

0.75 

-1214.18 

214.15 

4.71 

-1212.21 

212.21 4.71 

1.00 

-999.18 

0.00 

4.71 

-1000.00 

0.00 6.28 



Figure 3. Cantilever under end shear: exact and computed responses 


11.2 Cantilevered Beam Under Shear Load 


A cantilever beam is now subject to a vertical tip load. The results obtained with 16 2D 
Timoshenko TL beam elements in a term project by Abedzadeh, Mehrabi and Lofti 11 are 
compared in Figure 3 with the analytical solution given by Mathiasson et.al. s The 2D 
Timoshenko elements follow the exact solution without appreciable error. 


Remark 11.1. As noted in Remark 1.1, Mathiasson et. al. 8 reported fair to poor results beyond 
moderate rotations in this problem using a TL-based 2D Hermitian beam element. The difficulty 
can be traced to their use of an approximate expression for the curvature: 


'Uv' 

K % 

[l + (^) 2 f /2 ’ 

where primes denote derivatives with respect to X, instead of the correct small-strain 

K = — = o' — M v(l + ttx) - uWy 
dX [l+(«x) 2 + (^) 2 ] 3/2 ‘ 


( 11 . 1 ) 

TL value 
( 11 . 2 ) 
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Figure 4. Buckling of cantilever column: load-displacement responses 

The expression (11.2) rapidly losses accuracy as u' x and increase. Since (11.1) usually overes 
timates the actual curvature, it tends to overstiffen the element. 


11.3 Euler Buckling of Cantilever Column 

This buckling example is taken from Mathisen’s thesis. 12 The critical buckling load was 
traversed by treating the bifurcation point as a limit point by introducing a slight geometric 
imperfection. The following inputs were used in the study reported in Alexander et.al. . 
U = 100, A = 0.05, / = 0.20, E = 10 6 , G = E/2 and applied load P = 49.34A. The 
load response curves using the axial displacement as control variable are shown in Figure 
4. The computed results agree well with elastica solutions up to deflections of the order of 

the column length. 

11.4 Large Displacement of a 45° Cantilever Bend 

This 3D beam problem concerns the large displacement analysis of a 45° cantilever bend 
initially lying on the horizontal {X,Y} plane. The bend is an arc of a circle of radius 
r = 100 and the beam cross-section is a square with sides of unit length. The beam has 
modulus E = 10 7 and Poisson’s ratio v = 0. It is subjected to an end load P = 600 normal 
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x ( 29 . 3 , 70 . 7 , 0 . 0 ) 

Figure 5. Curved cantilever bend under tip load: Problem definition 
to the (AjT) plane as shown in Figure 5. 



Figure 6. Curved cantilever bend under tip load: Deflected shapes 

This problem was treated by Bathe and Bolourchi 13 with 3D brick elements, and sub- 
sequently by Simo and Vu-Quoc 14 and Cardona 15 with beam elements based on other 
formulations. Results with the TL 3D element described in Sections 5.5 and 10 were ob- 
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Table 4. Comparison of results for the 45° bend cantilever beam 


Source 

Load P = 300 

Load P = 450 

Load P = 600 

Bathe 13 

Simo 14 

Cardona 15 

Crivelli 5,7 

22.33, 58.84, 40.08 
22.50, 59.20, 39.50 
22.14, 58.64,40.35 
22.31, 58.85, 40.08 

18.62, 53.32,48.39 

18.38, 52.11, 48.59 
18.59, 53.34,48.39 

15.79, 47.23, 53.37 
15.90, 47.20, 53.40 
15.55, 47.04, 53.50 
15.75, 47.25, 53.37 



Length 

Polar Moment 
Second Moments 
Young’s Modulus 
Shear Modulus 


L = 240 mm 
I x = 2.16 mm 4 
ly - h = 0.0833 mm 4 
E = 71240 Nlmm2 
G = 27190 N/mm2 


Figure 7. Cable hockling: Problem definition 


tained by Crivelli 5 using 8 beam elements and applying the load in 6 equal increments. The 
solution method is incremental-iterative with full Newton iteration used in the corrective 
phase. Results for the three tip displacement components are compared with those of the 
aforementioned references in Table 4. Deflected shapes for selected load levels are shown 
in Figure 6. 

As can be observed the present results compare especially well with those obtained with 
3D elements in Bathe and Bolourchi 13 . 

11.5 Cable Hockling 

The second problem, cable hockling , is more challenging as regards modeling and post- 
buckling response analysis. An initially straight cable, modeled with 3D Timoshenko beam 
elements, is subjected to a tip torsional moment. The geometry and physical properties 
are given in Figure 7. The cable is clamped at one end and supported at the other so 


45 




that the only motions allowed at that point are axial displacement and torsional rotation 
about the longitudinal X axis. No rotation is allowed about Y or Z. The purpose of these 
restrictions is to keep the problem conservative, because if the torque is allowed to rotate 
about Y or Z, the problem becomes nonconservative and dynamical methods are required 
to assess its stability. 16,17 

This problem has received a great deal of attention from the engineering community due 
to its practical importance. The main objective is to estimate the critical applied torque at 
which the cable departs (bifurcates) from its straight configuration, resulting in the forma- 
tion of a loop or hockle. This has direct application to marine cables used in tasks such as 
lifting objects from the ocean floor, for which structural failure could be disastrous. Under 
the assumptions of infinitesimal bending deformations, Greenhill obtained an analytical 
formula to predict the critical torque; see e.g. pp. 417-418 of Love. 18 The post-bifurcation 
response analysis of this problem, however, has not been pursued until recently, as dis- 
cussed in the research conducted by Nour-Omid and Rankin. 19 This post-critical response 
has also been analyzed by using the present TL formulation. The structure is discretized 
by twenty equally- spaced beam elements. 

The deformed shapes at different load levels are shown in Figure 8, which displays the 
loop-formation process previously described. The curves on the left and right side show 
deformed shapes looking along the Y and Z axes, respectively. 

If the torque is held under the critical value, the beam twists without lateral deflection. 
Along this fundamental path the response is linear. At the critical torque a bifurcation 
point is reached, at which the fundamental path becomes unstable. The beam acquires 
a helical shape with the free end moving towards the clamped end. This new equilib- 
rium branch is unstable and the cable undergoes large displacements and rotations as the 
moment decreases. The unloading process continues until reaching a sharp limit point 
which corresponds to a negative value of the applied torque (it should be mentioned that 
Nour-Omid and Rankin, 19 who use a Hermitian beam element based on a corotational 
description, characterize this critical point as a secondary bifurcation; evidence for this 
classification is presently inconclusive). After traversing this point, the torque reverses 
again until the cable reaches a circular-shaped unloaded deformed configuration. 

The computed variation of the twist angle at the moving end versus the applied torque is 
given in Figure 9. Results are compared to those given by Nour-Omid and Rankin. 

12. CONCLUDING REMARKS 

This article covers an alternative technique for formulating geometrically nonlinear me- 
chanical finite elements based on the Total Lagrangian kinematic description. The Core- 
Congruential Formulation, or CCF, can be approached and studied at several levels of 
complexity. These give rise to what are called here the Direct, Algebraic-Generalized, and 
Differential- Generalized CCF. 

All of these variants, however, share a basic staged approach to the derivation of discrete 
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Figure 8. Cable hockling: Deformed shapes at different load levels: 

a) After bifurcation, b) M x = 50, c) M x = 0, d) M x = -50, 
e) At the limit point, f) M x = 0. 


finite element equations. In the innermost level, core equations are obtained at the particle 
level. These physically transparent equations depend only on the strain and stress measures 
adopted and the choice of terms retained in the strain energy. This is followed by a 
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Twist Angle, radians 

Figure 9. Cable hockling: Computed response as 
tip twist rotation vs. applied moment 

transformation phase that ultimately ends in the element nodal degrees of freedom. For 
complex elements the transformation phase is frequently carried out in stages. 

For elements whose nodal degrees of freedom include 3D finite rotations, a multistage 
transformation is convenient in the sense that the decision on which finite rotation mea- 
sure to use can be relegated to the next-to-last stage, while the choice of finite element 
interpolation and nodal degrees of freedom is introduced in the last stage. This strategy 
facilitate application of inner-level equations to sets of different but related elements and 
fosters programming modularity. 

At this point one may well pose general questions such as: Why use the Total Lagrangian 
description? Wouldn’t a co-rotational or Updated Lagrangian formulation be preferable? 
The answer is that each description has strengths and weaknesses. Here are advantages 
that can be cited for the TL description: 

1. If the element development can be carried out under the framework of the DCCF or 
AGCCF, a symmetric tangent stiffness formulation is guaranteed. This attribute has 
obvious advantages in stability analysis and traversal of bifurcation points. 

2. The choice of a fixed reference frame has advantages in nonlinear dynamic calcu- 
lations in that the mass matrix remains symmetric, with the same sparsity as in 
small-deflection analysis. 
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3. The use of Green- Lagrange strains and conjugate PK2 stresses, linked up with careful 
avoidance of hazardous kinematic approximations, automatically takes caxe of rigid 
body motions. No special filters to eliminate self-straining at the element level axe 
needed. 

4. Extension to nonlinear constitutive equations and finite strains is straightforward 
although may be laborious. In the case of the CCF, extension to small-strain material 
nonlinearity would affect only the core equations because the transformation phase is 
entirely governed by element kinematics. 

The TL description, however, suffers from several shortcomings: 

5. It makes no effective reuse of existing linear finite elements. This is the key strength 
of the element- independent corotational description. 19,20 

6. Has difficulties with the specification and computational treatment of boundary con- 
ditions intrinsically linked to the deformed configuration; for example pressure loads. 
Both the corotational and Updated Lagrangian descriptions handle this aspect better. 

7. Eventually breaks down for exceedingly large motions, for example rotations exceed- 
ing 2n. This disadvantage can be important for certain aerospace and mechanical 
structures. 

From this list it follows that there is no “kinematic description for all seasons.” One 
intriguing research area that may merit exploration in this regard is a combination of the 
TL and corotational descriptions that maintains their individual strengths while alleviating 
the more serious disadvantages. 
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Appendix 1 

Equivalence of DCCF and Standard TL Formulation 


The correspondence between the Direct Core Congruential Formulation (DCCF) and the Standard 
Formulation (SF) of the Total Lagrangian (TL) kinematic description is generally established for 
3D continuum finite elements. This connection was worked out in a course term project. 11 Such 
elements fit within the DCCF framework because their physical DOFs (node displacements) are 
of translational type. 

The Standard Formulation is based on the same scheme used for linear finite elements: first 
interpolate, then vary. As in the linear case, the departure point is extremization of the Total 
Potential Energy functional (TPE) over the element domain: 

J = U-W= f e T s 0 dV+\ f e T E edV - f u T b dV - f u T t dS, (Al.l) 

Jv 0 Jv, 0 Jv 0 Jsto 

where as usual conservative dead loading is assumed. In (Al.l), b is the prescribed body force field, 
t are surface tractions prescribed over portion Sto of the boundary in Co, and other quantities are 
as defined in Section 4. The weak equilibrium equations are obtained on making (Al.l) stationary: 

SJ = 6U-SW= f Se T s 0 dV + j 6e T EedV - j Su T bdV - f Su T tdS = 0. (A1.2) 

J Vq J Vq JvQ J S t Q 

The displacement and strain fields are interpolated in terms of the element degrees of freedom v: 

u = N v, 6u = N 8v, 6e = B£v, (A1.3) 

where B = B(v) depends in v but N does not. Substituting these interpolations into (A1.2) 
yields the residual equilibrium equations 

6 J = 6v T r = 6v T (f - p) = 0. (A1.4) 


where 


f = 



B r (s 0 +E e)dV = 


/ B T sdV', 

p = j N r b dV+ I Ns dS, 

(A1.5) 

'Vo 

Jv o Js t o 



where f and p are the internal and external force vectors, respectively, and s = s° + Ee are the 
PK2 stresses in C. Because the variations 8v are arbitrary, the residual-force nonlinear equilibrium 
equation isr = f— p = 0orf = p. The tangent stiffness matrix is given by 


K = 


dr _ df 
dv dv' 


(A1.6) 


because p (for conservative dead loading) does not depend on v. Splitting B — B c + B„(v), where 
B c is constant but B v depends on v, gives the well known decomposition 


K = K 0 + K d + K g , 


(A1.7) 
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where Ko, K d and Kg denote the linear, initial-displacement and geometric stiffness matrices, 
respectively. These are given by 


Ko = f 

Jv Q 


BjEBcdV, 


K D = I (Bf EB„ + B^EB C + B^EB„) dV, 
Jv 0 

fv 0 


(A1.8) 


K g Sv 


6B t s dV. 


To correlate these standard forms with those produced by the DCCF, we note that the GL strains 
can be also split as e = e c + e v , where e c and e„ are linear and nonlinear in v, respectively. The 
latter may be expressed in terms of the displacement gradients as 


e„ = ±Ag, 


(A1.9) 


where A is the 6x9 matrix 


A = 


■gf 

0 

0 - 


'9 1 

92 

9 3 

0 

0 

0 

0 

0 

0 - 

0 

si 

0 


0 

0 

0 

94 

9 5 

9 6 

0 

0 

0 

0 

0 

g-T 


0 

0 

0 

0 

0 

0 

9 7 

9s 

9 9 

0 

gj 

si 


0 

0 

0 

9 7 

9s 

59 

94 

9 5 

96 

gj 

0 

si 


9 7 

9s 

9 9 

0 

0 

0 

9 1 

92 

9 3 

Lg? 

sT 

0 . 


-94 

9 5 

9 6 

9 1 

9 2 

9 3 

0 

0 

0 . 


(A1.10) 


in which the displacement gradients are vector-arranged as 

g r = [<7i 92 ••• gs ixt ■“ 


Comparing 


6e v = ^ 6 A g + i A <5g = A 6g, 


(Al.ll) 


(A1.12) 


to the DCCF transformation relation <5g = G(5v, in which G is independent of v, we see that 

B, = AG. (A1.13) 

The other expression we require is 6 A T s, which appears in the geometric stiffness matrix con- 
tracted with 6v: 


G t 6A T sdV. 


K g Sv= / SB T sdV= I 
Jv 0 JVq 

It is well known — see for instance Chapter 19 of Zienkiewicz 21 — that 


(A1.14) 


SA. s = M Sg — MGSv, with M = 


Sl I S4I S5I 

S 4 I S 2 I Sgl 

LssI S6 1 s 3 l 


(A1.15) 
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where I is the 3x3 identity matrix and s,-, i = 1,...6 are components of the PK2 stress tensor 
ordered = sn, S 2 = S 22 , ... 56 = $ 23 . Using this relation, Kg can be placed in the standard 
form 


K g 


•fv 0 


G t M GdV, 


(A1.16) 


which by inspection is seen to be the DCCF-transformation of the core geometric stiffness M = 
S g = 5, Hi, with the H, matrices defined in (4.3). 

To correlate other terms, write the linear part of the GL strains in terms of gradients as 


e c = Dg = DG 6v, with D = 


-1 

0 

0 

0 

0 

0 

0 

0 

0 - 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

1 

0 

1 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

1 

0 

0 

.0 

0 

0 

0 

0 

1 

0 

1 

0 . 


(A1.17) 


The numerical D matrix can be easily related to the h< vectors introduced in (4.3). Because 
both D and G are independent of v it follows that 6e c — DG£v and consequently B c = DG. 
Partitioning A as [af ... a<f ] one easily finds that a, = Hig. Now the following identities 

can be verified through simple algebra: 


D t ED = Eijhihf = So, 

D t EA = En hmj = Eij h,g T H; = Sj , A t ED = Sf, 

A t EA = Eija< a ; = ^ijHigjgfHj = S 2 = S^, 

M = s- Hi + Ei/higH; + |(g T Hig)Hj = s?H< + S* + = s.H,. 


(A1.18) 


Comparing these to the expressions of Section 4.3 we conclude that 



(A1.19) 


which displays the equivalence of both formulations when no approximations are made. This proof 
may be extended without difficulty to the AGCCF in which case G is a function of v, although 
as noted in the text that situation is sometimes mishandled in the Standard Formulation through 
the introduction of a priori kinematic approximations. The equivalence between DGCCF and SF 
is more difficult to prove because there is no TPE functional from which the latter can be derived, 
and such connection should be regarded as an open problem. 
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Appendix 2 

Tangent Stiffness Subroutine for Two-Node 3D Bar 


The subroutine listed below implements the expression derived in Section 7.1 for the global tan- 
gent stiffness of a two-noded TL 3D bar element. The Fortran implementation is biased in favor 
of computational speed. Therefore, it contain no loops or calls to other subroutines, as that 
would slow down the calculations. The formation of the 6x6 tangent stiffness matrix requires 74 
multiplications, 14 divisions, 64 additions/subtractions and 3 square roots, for a total of approx- 
imately 160 double precision floating-point operations. On a 15-MFlop workstation, this results 
in approximately 50000 elements formed per CPU second. 

One point that deserves some attention is the numerically-stable choice of the local coordinate 
system. Axis X [s uniquely defined as the longitudinal bar axis that passes through the end nodes, 
but axes Y and Z may be arbitrarily rotated about X because a 3D bar has no preferred transverse 
directions. Although the end result, namely K, is independent of this choice it is important to 
choose Y and Z in a stable manner that works for any element orientation. The procedure is 
described in code comments below. 


C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 


subroutine BAR3K (xO, v, em, area, sO, k, status) 
Compute tangent stiffness matrix of 2-node TL 3D bar 


Inputs: 

xO 

v 

em 

area 

sO 


Reference node coordinates 
Global node displacements 
Elastic modulus 

Cross section area in reference configuration 
PK stress in reference configuration 


Outputs : 
k 

status 


Tangent stiffness of bar element 

Blank if no error detected else error message 


double precision 
character* (*) 
double precision 
double precision 
double precision 
double precision 
double precision 
double precision 
equivalence 
equivalence 
equivalence 


x0(3,2), v(3,2), em, area, sO, k(6,6) 
status 

xbar(3) , ybar(3) , zbar(3) 
lxbar, lybar, lzbar, lxbar2, 1x2 
dx(3) , du(3) , g(3), tt(3,3), ea, sg, s 
ktll, ktl2, kt 13 , kt21 , kt22, kt23 
kt31 , kt32 , kt33 

sbarll , sbarl2 , sbarl3, sbar22, sbar23 , sbar33 
(tt(l,l) ,xbar(l)) 

(tt (1 ,2) ,ybar(l) ) 

(tt(l ,3) ,zbar(l)) 


status = 


Form 3x3 global to local transformation matrix T 
Note: array tt receives T-transpose (inverse of T) 

Method: begin by getting dir cosines tll,tl2,tl3 of Xbar vrt X 
Compute sqrt(tll**2*tl2**2); if gt tol set directors of Ybar 
to -t 12 , t 11 ,0 and normalize; else set to 0,-tl3,tl2 and normalize. 
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Value of tol set to .01 but any value ge 0 would do fairly well. 
Finally Zbar = Xbar x Ybar . 

dxCl) = x0(l ,2) - x0(l , 1) 
dx(2) = x0(2 f 2) - x0(2 , 1) 

dx(3) = x0(3 ,2) - x0(3 , 1) 

lxbar2 = dx(l)**2 + dx(2)**2 + dx(3)**2 

lxbar = sqrt(lxbar2) 
if (lxbar .eq. 0.0) then 

status = ’ BAR3K : Bar has zero length * 

return 
end if 

xbar(l) = dx( l)/lxbar 

xbar(2) = dx(2)/lxbar 

xbar(3) = dx(3)/lxbar 

lybar = sqrt (xbar (1) **2 + xbar(2)**2) 
if (lybar .gt. 0.01) then 

ybar(l) = -xbar (2) /lybar 
ybar(2) = xbar ( 1) /lybar 

ybar(3) = 0.0 

else 

lybar = sqrt (xbar(2)**2 + xbar(3) **2) 
ybar(2) = -xbar (3)/lybar 
ybar(3) = xbar (2) /lybar 

ybar(l) = 0.0 

end if 

zbar (1) = xbar (2) ♦ybar (3) - xbar (3) *ybar(2) 
zbar (2) = xbar (3)*ybar(l) - xbar(l)*ybar(3) 
zbar (3) = xbar (1) *ybar(2) - xbar(2)*ybar(l) 
lzbar = sqrt (zbar ( 1) **2 + zbar(2)**2 + zbar(3)**2) 
if (lzbar .ne. 1.0) then 

status = y BAR3K : lzbar ne 1.0: cannot happen’ 
return 
end if 

Form core material and geometric stiffness in local coordinates 
Note: only six different entries of sbax need to be computed 
because 

R -R 

C Kbar= 


c 


-R R 




c 

where 





c 


sbar 1 1 

sbarl2 

sbarl3 


c 


R = 

sbar22 

sbar23 


c 

n 


symm 


sbar33 



ea = 

em*area/lxbar 






du(l) = 

v(l ,2) - v(l,l) 





du(2) = 

v(2 ,2) - v(2 , 1) 





du(3) = 

v(3 ,2) - v(3 , 1) 





1x2 = 

(dx(l)+du(l) )**2 + (dx(2)+du(2))**2 + (dx(3)+du(3) ) **2 



sO + em*0.5*(lx2-lxbar2)/lxbar2 


sg = s*area/lxbar 

g ( 1 ) * (tt(l, l)*du(l)+tt(2, D*du(2)+ttC3, l)*du(3)) /lxbar 
g(2) = (tt (1 ,2) *du(l) +tt (2,2) *du(2)+tt (3 ,2)*du(3) )/lxbar 
g(3) = (tt (1 ,3)*du(l)+tt(2 f 3)*du(2)+tt(3,3)*du(3))/lxbar 
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sbarll = ea * 
sbar22 = ea * 
sbar33 = ea * 
sbarl2 = ea * 
sbarl3 = ea * 
sbar23 = ea * 


(1 .+g(l))**2 + sg 
g(2)**2 + sg 

g(3)**2 + sg 

(l.+g(l))*g(2) 
(l.+g(l))*g(3) 
g(2)*g(3) 


Transform 

core stiffness to physical global coordinates 

ktll = 

tt(l , l)*sbarll + tt(l ,2)*sbarl2 + tt(l,3)*sbarl3 

kt21 = 

tt (1 , l)*sbarl2 + tt(l,2)*sbar22 + tt(l ,3)*sbar23 

kt31 = 

tt(l, l)*sbarl3 + tt(l,2)*sbar23 + tt(l ,3)*sbar33 

ktl2 = 

tt (2, l)*sbar 11 + tt(2,2)*sbarl2 + tt(2,3)*sbarl3 

kt22 = 

tt(2, l)*sbarl2 + tt(2,2)*sbar22 + tt(2,3)*sbar23 

kt32 = 

tt(2, l)*sbarl3 + tt(2 ,2)*sbar23 + tt(2,3)*sbar33 

kt 13 = 

tt(3,l)*sbarll + tt(3,2)*sbarl2 + tt(3,3)*sbarl3 

kt23 = 

tt(3, l)*sbarl2 + tt(3,2)*sbar22 + tt(3,3)*sbar23 

kt33 = 

tt(3, l)*sbarl3 + tt(3,2)*sbar23 + tt(3,3)*sbar33 

k(l , 1) = 

tt(l,l)*ktll + tt(l,2)*kt21 + tt(l,3)*kt31 

k(l,2) = 

tt(l,l)*ktl2 + tt(l ,2)*kt22 + tt(l,3)*kt32 

k(2 , 1) = 

k(l,2) 

k( 1 ,3) = 

tt(l , l)*ktl3 + tt(l ,2)*kt23 + tt(l,3)*kt33 

k(3 , 1) = 

k(l ,3) 

k(2 ,2) = 

tt(2, l)*ktl2 + tt (2 , 2) *kt22 + tt(2,3)*kt32 

k(2,3) = 

tt(2,l)*ktl3 + tt(2,2)*kt23 + tt(2,3)*kt33 

k(3 ,2) = 

k(2 ,3) 

k(3,3) = 

tt(3 , l)*ktl3 + tt(3,2)*kt23 + tt(3,3)*kt33 

k(4, 1) = 

-k(l , 1) 

k(5, 1) = 

-k(2 , 1) 

k(6 , 1) = 

-k(3 , 1) 

k(4,2) = 

-k(l ,2) 

k(5,2) = 

-k(2,2) 

k(6 , 2) = 

-k(3 ,2) 

k(4,3) = 

-k(l ,3) 

k(5,3) = 

-k(2,3) 

k(6,3) = 

-k(3,3) 

k(l,4) = 

k(4, 1) 

k(2 ,4) = 

k(4,2) 

k(3 ,4) = 

k(4,3) 

k(4 ,4) = 

kCl.l) 

k(5,4) = 

k(2 , 1) 

k(6,4) = 

k(3, 1) 

k(l ,5) = 

k(5, 1) 

k(2,5) = 

k(5 ,2) 

k(3 ,5) = 

k(5 ,3) 

k(4,5) = 

k(5 ,4) 

k(5,5) = 

k(2 ,2) 

k(6 ,5) = 

k(3,2) 

k(l ,6) = 

k(6, 1) 

k(2,6) = 

k(6,2) 

k(3 ,6) = 

k(6 ,3) 

k(4,6) = 

k(6 ,4) 

k(5 ,6) = 

k(6 ,5) 

k(6 ,6) = 

k(3,3) 

return 


end 
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